
Calhoun 

iniQiuiic^iul Ar{hiv« of tilt Mil vdl Poii^roduiit School 


Calhoun: The NPS Institutional Archive 
□Space Repository 



Theses and Dissertations 


1. Thesis and Dissertation Collection, all items 


1993-09 

Comparison of higher order moment spectrum 
estimation techniques 

McAloon, Jeffrey F. 

Monterey, California. Naval Postgraduate School 


http://hdl.handle.net/10945/39976 


This publication is a work of the U.S. Government as defined in Title 17, United 
States Code, Section 101. Copyright protection is not available for this work in the 
United States. 

Downloaded from NPS Archive: Calhoun 



DUDLEY 

KNOX 

LIBRARY 


htt p://w ww. n ps. e du/l ib ra ry 


Caflwuo is the Naval Postgraduate School's public access digital repository for 
research mate rials and institutiional publicatkins created by the NPS community. 
Calhoun is named for Professor of Mathematics Guy K. Caftiouo, NPS's first 
appointed — and published — schoteily author. 

Dudley Knox Library / Naval Postgraduate School 
411 Dyer Road / 1 Univefsity Circle 
Monterey, California USA 93943 






UJNCLASSlMiiU 

SECUf=tlTY CLASSIFICATION OF THIS PAGE 




REPORT DOCUMENTATION PAGE 


UBITY CLASSIFICATION UNCLASSIFIED I RESTRICTIV 


2b. DECLASSIFlCATIONA)OWNGRADING SCHEDULE 


4. PERFORMING ORGANIZATION REPORT NUMBER(S) 


3. DISTRIBUTION/AVAIUBILITY OF REPORT 
Approved for public release; 
distribution is unlimited 


5. MONITORING ORGANIZATION REPORT NUMBER(S) 



6c. ADDRESS (CJIy, State, and ZIP Code) 

Monterey, CA 93943-5000 




7a. NAME OF MONITOBI 

Naval Postgraduate School 


7b. ADDRESS (City, State, and ZIP Code) 

Monterey, CA 93943-5000 



8c. ADDRESS (City, State, and ZIP Code) 


9. PROCUREMENT INSTRUMENT IDENTIFICATION NUMBER 


10. SOU 




WORK UNIT 
ACCESSION NO, 


11. TITLE (Include Security Classification) 

COMPARISON OF HIGHER ORDER MOMENT SPECTRUM ESTIMATION TECHNIQUES (U) 


AL AUTHOR{S) 


McAloon, Jeffrey F. 




14. DATE OF REPORT (Year, Month, Day) 


esis 

FROM TO. 

September 1993 


The views expressed in this thesis are those of the author and do not reflect the official 
policy or position of the Department of Defense or the United States Government. 


18. SUBJECT TERMS (Continue on reverse if necessary and identify by block number) 

bispectrum; 1-1/2 D; IPS; higher-order moments; cumulants 


17. 

COSATI CODES 

FIELD 

GROUP 

SUB-GROUP 1 


19. ABSTRACT (Continue on reverse if necess^ and identify by block number) 

This thesis compares the detection performance of the 1-1/2 D instantaneous power spectrum (1-1/2 Djpj), the 
bispectrum, the instantaneous higher-order moment slice (IHOMS) method, and the spectrogram for 
multi-component stationary signals, harmonically related stationary signals, and multi-component linear FM signals 
corrupted by additive white Gaussian noise. In addition, a determination of the relative processing gain between the 
1-1/2 Dips method and the spectrogram is made for stationary signals in noise. 

The results of this thesis show that 1-1/2 Dips ® processing gain advantage over that of the spectrogram for a 
range of input SNR that depends upon the size of the data window. Under some conditions, the bispectrum can 
detect both harmonic coupling and phase coupling between the components of multi-component signals. IHOMS’ 
ability to detect linear chirps in noise is limited to chirps having different slew rates, and the method has a significantly 
greater computational cost than both the spectrogram and 1-1/2 Djpg. 


20. DISTRIBUTION/AVAILABILITY OF ABSTRACT 21. ABSTRACT SECURITY CLASSIFICATION 

g unclassified/unlimited Qsameasrpt. Qdtic users UNCLASSIFIED 




OD FORM 1473,84 MAR 83 APR edition may be used untH exhausted SECURITY CLASSIFK 

All other editions are obsolete 


SECURITY CLASSIFICATION OF THIS PAGE 

UNCLASSIFIED 




































Approved for public release; distribution is unlimited. 


COMPARISON OF HIGHER ORDER MOMENT 
SPECTRUM ESTIMATION TECHNIQUES 

by 

Jeffrey F. McAloon 
Lieutenant, United States Navy 
B.S.E.E., University of South Florida, 1983 


Submitted in partial fulfillment of the 
requirements for the degree of 


MASTER OF SCIENCE IN ELECTRICAL ENGINEERING 


from the 


NAVAL POSTGRADUATE SCHOOL 
September 1993 




u 




ABSTRACT 


This thesis compares the detection performance of the 1-1/2 D instantaneous power 
spectrum (1-1/2 Dips), bispectrura, the instantaneous higher-order moment slice 
(MOMS) method, and the spectrogram for multi-component stationary signals, 
harmonically related stationary signals, and multi-component linear FM signals corrupted 
by additive white Gaussian noise. In addition, a determination of the relative processing 
gain between the 1-1/2 Dips ni®thod and the spectrogram is made for stationary signals in 
noise. 

The results of this thesis show that 1-1/2 Dips ® processing gain advantage over 
that of the spectrogram for a range of input SNR that depends upon the size of the data 
window. Under some conditions, the bispectrum can detect both harmonic coupling and 
phase coupling between the components of multi-component signals. MOMS’ ability to 
detect linear chirps in noise is limited to chirps having different slew rates, and the method 
has a significantly greater computational cost than both the spectrogram and 1-1/2 Djps. 
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I. INTRODUCTION 


A. OVERVIEW 

Spectral estimation is an important data analysis tool used in signal processing to 
determine the distribution of power as a function of frequency. The “classical” Fourier 
transform methods and a variety of more recently developed parametric methods are now 
commonly used to estimate power spectral density. These methods are based on the 
autocorrelation domain and therefore rely upon information contained in the second order 
moment of the data sequence. Recent improvements in computational capability have 
opened up a new and active research area concerned with the extraction of additional 
information contained in the data sequence’s higher-than-second order statistics. [Refs. 
1 . 2 ] 

Higher order statistics involve cumulants rather than moments. Cumulants and 
moments are similar, and each can be expressed in terms of the other. In fact, under 
certain conditions they are identical. Higher order spectra (HOS), also known as 
polyspectra, are formed by taking the Fourier transform of the cumulant sequence. For 
example, the two-dimensional Fourier transform of the third order cumulant yields the 
polyspectrum known as the bispectrum. Likewise, the three-dimensional Fourier 
transform of the fourth order cumulant is called the trispectrum. If the first moment is 
zero, the second order cumulant is equal to the second order moment, and the bne- 
dimensional Fourier transform of either sequence produces the familiar power spectrum. 
[Refs. 1,3] 

Higher order statistics and their spectra possess some exploitable properties that their 
second order counterparts do not. One potential advantage involves a greater resistance 
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to the effects of certain types of additive noise. For white or colored, zero mean Gaussian 
processes, the cumulant of any order greater than two is zero and the corresponding 
polyspectrum is zero. Consequently, polyspectra are expected to be more resistant than 
the power spectrum is to the effects of additive Gaussian noise. One limitation of the 
power spectrum is that phase cannot be accurately reconstructed unless the underlying 
signal or system is minimum phase. In contrast, higher order spectra preserve true phase 
information. Finally, HOS are useful for detection and classification of non- 
linearities. Depending upon the application, the extra computational cost of HOS may be 
justified. [Refs. 3,4] 

The bispectrum, the 1-1/2 D instantaneous power spectrum [Ref. 5], and the 
instantaneous higher order moment sKce (IHOMS) method [Ref. 6] are studied in this 
thesis. Performance comparisons are made with respect to each method as well as to an 
accepted second order method. 

B. THESIS OUTLINE 

The following describes the organization of the remainder of this thesis. Chapter n 
first defines cumulants and shows their relation to moments. The essentials are then 
presented for each spectral estimation method studied. Chapter HI makes a processing 
gain comparison between the 1-1/2 D instantaneous power spectrum and the spectrogram. 
Chapter IV displays simulation results that show how the different methods perform on a 
variety of signals. Chapter V presents conclusions and suggestions for future 
work. Appendices contain the detailed steps of Chapter III calculations and the Matlab 
programs used to execute the computer simulations. 
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II. ELEMENTS OF HIGHER ORDER SPECTRAL ANLYSIS 

This chapter describes the higher order based techniques compared later in this 
thesis. Sufficient information is presented in order to develop an adequate understanding 
of each method. We first define cumulants and show how they are used instead of 
moments to form polyspectra. 

A. CUMULANTS AND MOMENTS 

In general, the characteristic function, O ( 7 ©). is defined as the conjugated Fourier 
transform of a random variable’s probability density function, f(x) [Ref. 7]: 

00 

d)(y©) = (f(x)exp(J(i}x)dx. (2.1) 

wOO 

The right-hand side of (2.1) is simply the expectation of exp (JdiX). The n^ moment of 

the random variable X can be found by evaluating the n^ derivative of the characteristic 
function with © set equal to zero [Ref. 7]: 

E[X"} = . (2.2) 

„.o 

In comparison, the n**^ cumulant is defined as the nth derivative of the natural logarithm of 
the characteristic function evaluated at © = 0 [Ref. 4]; 

C{X"} = (-;)«-^n[a>(©)] . (2.3) 

dor n 

Although (2.3) clearly shows the basic difference between a moment and a cumulant, 
it does not readily show that a cumulant can be considered as a moment with its lower order 
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statistical dependence removed [Ref. 8], This is perhaps best seen by computing 
cumulants in an alternate manner that involves a partitioning of equal and lower order 
moments. 

Given a set of random variables X - {Xp^ 2 , , Ix - {1,2 .is 

defined to be the set of component indices contained in X. Denoting a subset of Ix by I, 
using TTixU) to represent the moment of those components of X comprising Xj, and 
denoting the cumulant of Xj as Cx U ). the moment-to-cumulant (M-C) equation can be 
written as [Ref. 3]: 

where Ip are non-intersecting non-empty subsets of I that form the partitions, and q is var¬ 
ied fi-om one to the number of elements in I. The summation in (2.4) is over all unique 
partitions and for aU q. The first order cumulant is equal to the first order moment since 
only one unique partition exists for one random variable. Table 2.1 shows how (2.4) is 
applied to compute the moment partitions that form the second order cumulant. Table 
2.2 displays the details of the third order cumulant calculation. Summing the last column 
of Table 2,1 yields the second order cumulant. Summing the last column of Table 2.2 
produces the third order cumulant. The first three orders of cumulants, expressed in 
terras of the component indices of a set of random variables, are then: 
c(l) = m(l), 

c(l,2) = m(l,2)-m(l)m(2), 

c(l,2,3) = m(l,2,3)+2m(l)m(2)m(3) ~m(l,2)m(3) 
-/7j(l,3)m(2 )-m(2 ,3 )/7 z(l). 

Calculation of the fourth order cumulant is similarly detailed in [Ref. 3]. An expression 
that calculates moments from partitioned cumulants is also given in [Ref. 3]. 
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TABLE 2.1: COMPUTATION OF SECOND ORDER CUMULANT 


h 

h 

q 

M-C 

equation 

(2.4) 

1.2 


1 


1 

2 

2 

-m(l)m(2) 


TABLE 2.2: COMPUTATION OF THIRD ORDER CUMULANT 


h 

h 

h 

q 

M-C equation (2.4) 

1,2.3 



1 

m(l,2,3) 

1.2 

3 


2 

-m(l,2)m(3) 

1.3 

2 


2 

-m(l,3)m(2) 

2.3 

1 


2 

-m(2,3)m(l) 

1 

2 

3 

3 

2m(l)m(2)m(3) 


Inspection of the moment-cumulant expressions show that dependencies among 
random variables are removed when cumulants are calculated. In fact, both the second 
and third order cumulants are zero if all the random variables are independent of each other. 
[Ref. 8] 

Another observation regarding the above moment-cumulant expressions applies in the 
commonly encountered situation where the means of the random variables are equal to 
zero. When this is true, the first order cumulant is zero since it equals the mean, the 
second order cumulant simply equals the variance, and the third order cumulant equals the 
third order moment. [Ref. 3] 
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There is more than one acceptable choice when it comes to choosing which term(s) to 
conjugate in the expectation expression for moments of order greater than two. Different 
choices have different consequences when the signal being processed is complex [Refs. 9, 
10]. One consequence is considered later when the symmetry region of the bispectral 
plane is discussed. For now, the conjugated terms are chosen arbitrarily. The following 
expressions show the details of the moment calculations in the last set of equations for the 
zero mean case [Ref. 11]: 

= E{x[n]} =0, (2.5) 

=£{/[n]x[n + /i]}, (2.6) 

[/i./j] = £{/[n]x[/j + /ilx[n + /2] } . (2.7) 

The expression for the fourth order cumulant as calculated by (2.4) has 15 moment terms. 
If the mean is zero all but four terms go to zero when the signal is real [Refs. 3,11]: 

[/j, I2, ^3] = £■ {x [rt] x[n + /j] X [n +12] x [n + /j] } 

~E {x [n] X [/i + /^]} £ {x [rt + 12 ] X [n + 13 ] } 

-£■ {x [n] X [aj +/ 2 ] } £■ {x [n +/j]x [n +/g] } (2.8) 

-£ (x [n] X [n + 73] } £■ {x [n + /j] X [n + /2] } . 

Provided that two of the four random variables are conjugated in the expectation terms 

of (2.8), the fourth order cumulant of a zero mean complex process has just three 

expectation terms. One of the last three terms in (2.8) will be zero. The zero term 

depends upon which variables arc conjugated. [Ref. 11] 

Cumulants have three properties that make them more desirable than moments when 

it comes to higher order statistics: 

1. Each cumulant is independent of all lower order cumulants. Consequently, all 
cumulants of order greater than two are equal to zero for a Gaussian process, as 
a Gaussian process is completely characterized by its first and second 
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moments. Higher order moments, on the other hand, can contain information 
about lower order moments. In fact, only higher order moments of odd order 
are zero for a Gaussian process. The non-zero even order higher moments do 
not contain any information not already contained in the first two moments. 
[Ref. 8] 

2. If a set of n random variables can be divided into more than one statistically 
independent subset, then the n* order cumulant, unlike the n**^ order moment, is 
equal to zero. [Ref. 4] 

3. Unlike moments, the cumulant of the sum of two independent stationary random 
processes is equal to the sum of the cumulants of each process. [Ref. 4] 

B. THE BISPECTRUM 

1. Definition and Region of Support 

The bispectrum is defined as the two-dimensional Fourier transform of the third 
order cumulant [Ref. 4]: 

00 oo 

S(«>l-®2) =7^2 E £ c‘”(i„/2)«p{-;(0),;i+m2y}. (2.9) 

Analogous to the power spectrum, the bispectrum can also be defined with frequency 
domain quantities. Given N samples of stationary signal x (n), its Fourier transform is 

N-l 

X(©) = ^x(n)exp(-j(on) . (2.10) 

n-0 

Assuming that the third order cumulant in (2.9) is computed with the conjugation scheme 
shown in (2.7), the equivalent frequency domain expression for the bispectrum is obtained 
through an extension of the periodogram [Refs. 9,12]; 
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(2.11) 


fi (co^, (Dj) * (oo^) X (t02) X* (<0j + 002) ; 

where |co^| ^ n, |(D 2 | ^ n, and jto^ + ( 02 j ^ n. Equations (2.9) and (2.11) represent two 

different non-parametric methods that can be used to calculate the bispectrum. The 
approach taken by (2.9) is called the indirect method, while the approach used by (2.11) is 
known as the direct method. [Ref. 1] 

In general, the bispectrum’s region of support is a hexagon centered at the origin 
of the (cOj, £ 02 ) bi-frequency plane. Evaluation of the mean product of three Fourier 

amplitude terms shows that the bispectral plane exhibits certain symmetries. For a 
stationary, real, continuous time signal the expected value of the product of three Fourier 
components is [Ref. 13]: 

£{X^(©j)X^(a)2)Xc(t03)} = B^((D^,©2)5 (o>i + “2‘''®3)* 

where the subscript c denotes continuous time quantities. 

Three properties of (2.12) determine bispectral symmetry. First, since the 
frequency indexes must sum to zero, ©3 = - ©^ *- © 2 . Second, the frequency indexes in 
the expectation operation can be interchanged. Third, for the bispectrum of a real signal, 
conjugate symmetry results since B(-©j,-© 2 ) =^^(©^,© 2 ). Using the first 
property to express ©3 in terms of ©^ and © 2 , the symmetry lines for a continuous time 
signal are [Ref. 13]; 

©j » ©2, 

2©j = -©2 (from ©j 5s ©3), 

2©2 = “©1 (from ©2 = © 3 ), 

©1 = -© 2 , 

©2 = 0 (from ©3 = -© 3 ), 
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CO, = 0 


(from GO2 - “0>3)- 


271 


If the signal is band-limited to and sampled at a rate of = 2©^ = , the 

. 2nn Inn , ■ 

bispectrum term in (2.12) becomes B(©j + —-.© 2 +after applying the 


approximation; 


_ ^ 2nn^ 

X(©)s £ x,(© + —). 

nm~eo 

Inn 


(2.13) 


The frequency indexes in (2,12) now must sum to instead of zero. As a result, the 

following set of symmetry lines are formed in the bispectrum of a sampled signal [Ref. 
13]; 


©j = © 2 . 


271/1 

2©j +©2 = 


2nn 

©^+2©2 = 


©^ = “®2’ 


2nn 

©2 = 


2nn 

m, s= - 

1 T 


Figure 2.1 shows the region of support and its symmetry lines for the first of an 
infinite number of bispectral hexagons given a sampled real signal [Refs. 12,13]. The 
sampling period is assumed to be equal to one. A complex signal has the same hexagonal 
support region as that of a real signal but it is symmetrical about only one of three dotted 
symmetry lines shown in Figure 2.1. Depending upon the term conjugated in (2.7) used to 
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compute the third order cumulant for the indirect method, or which term in (2.11) is 
conjugated to compute the bispectrum by the direct method, the bispectrum of the complex 
signal exhibits symmetry in only one quadrant about either the co^ = rOj 

0)2 “ -2o)j line, or the = “2©2 line. The conjugation scheme used exclusively in 
this study is shown in (2.7) and (2.11). Under this scheme, the bispectrum of a complex 
signal displays two-fold symmetry in the first quadrant about the = ©2 line. There 
are two other possible ways to conjugate just one term, and three possible ways to conjugate 
any two terms in the applicable expressions. The bispectrum symmetry for these other 
schemes are considered fully in [Ref. 9]. 
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Because of the symmetry described above, only a small portion of the bispectrum 
needs to be computed. Once the bispectrum is known for this small region known as the 
non-redundant region, the rest of the bispectrum can be found using symmetry. The non- 
redundant region, again for the sampled signal case, is shown in Figure 2,2. [Ref. 12,13] 



Figure 2.2: The non-redundant region of the bispectrum. 

2. Computation by the Indirect Method 

The indirect method requires a two-dimensional Fourier transform of an 
estimated third order cumulant sequence. A two-dimensional window can be used 
provided that it satisfies certain requirements based on cumulant symmetries. The 
following procedure outlines the indirect method given the data set 
{A'(0).X(l)....A:(lV-l)};[Ref.l] 

1. Form K segments of M samples such that N=KM. 

2. Remove the mean from each segment. 

3. Defining i=l,2,...K to indicate the segment number, estimate the third order 

moment sequence for each segment, {jt ^ = 0 ,1,...M - 1}: 
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(2.14) 


1 * • • 
r^'\m,n) = x''^hl + m)x^'hl + n)-, 

where sl=max{0,-ni,-n) and s2=min(M-l,M-l-m, M-l-n). This sequence is 
also called the tri-correlation sequence. 

4. Average all segment moment estimates: 

Ri.m,ri) = (2.15) 

5. Multiply the averaged moment estimate by a two-dimensional window, W(tn,n), 
and take the two-dimensional Fourier transform to obtain the bispectrum esti¬ 
mate: 

L L 

B(c0j,a>2) = V y R{m,n)Wim,n)exp{-j{(£)^m + ai2n)}; (2.16) 

m^-L nm-L 

where L<M-\. 

The window, W{m,n), used in the tri-correlation sequence defined in (2.16) must 
satisfy the following conditions [Ref. 1,11]: 

1. The window conforms to third order cumulant symmetry. Specifically, 

W(m,n) = W{n,m) = W(-m,n-m) = W {m - n,-n) . 

2. The window is zero outside the region of support for the tri-correlation 
sequence. 

3. The window is normalized to one at msm=0. 

4. The Fourier transform of the window is non-negative. 

Windows used in this study belong to a separable class of window functions defined by 

W{m,n) ^ d{m)d(n)d{n-m), (2.17) 

where an appropriate one-dimensional window, d(m), was chosen subject to the following 
constraints: 
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d(m)=d(-m) 
d(m)=0 for m>L, 
d(0)=l, 

D(g))^ 0 for all CO. [Ref. 11] 

The unit hexagonal window was constructed by applying (2.17) to a one- 
dimensional rectangular window [Ref. 14], 

d{m) =1 for|m|£L, 
d{m) = 0 otherwise. 

The one-dimensional Parzen window, given by 

2 3 

d(m) * l-6(™) for|m|^|. 

rf(m) =2(1-^) for^^ImKL, 
d(m) s= 0 for|m|>L, 

was similarly transformed into a two-dimensional Parzen window [Ref. 1]. The last 
separable window considered is known as the minimum bias supremum window, or 
simply the optimum window. The corresponding one-dimensional window is defined by 
[Ref. 1]: 


d(m) = ^ 
7t 

d(m) = 0 


. .nm 
sin (—) 


ImL 

+ (l-^)cos(—) 


for \m\>L. 


for \m\ £L, 


A trade-off situation exists between the Parzen window and the optimum window 
with regard to the windows’ statistical properties. The bias of the optimum window is 
about 18% smaller than that of the Parzen window. However, the optimum window has 
a 26% greater variance. [Ref. 15] 
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3. Computation by the Direct Method 

The direct method generates the bispectrum from (2,11). This requires a one¬ 
dimensional Fourier transform of the signal. Frequency domain smoothing can be used 
to reduce the variance of the final bispectrum estimate. The direct method is outlined for 
thedataset {X(0),X(1), 1)} [Ref. 1,12]: 


1. The data sequence is segmented into K segments. The mean is removed from 
each segment before zero padding the segment to a convenient fast Fourier 
transform length, L. 

2. An L-point Fourier transform, (A.), is computed for each segment. The 
superscript designates the segment number. 

3. The bispectrum estimate of the segment, denoted as 5,- is then computed by: 


4. Two different approaches can be taken now to reduce the variance of the final 
estimate. The first approach, denoted Bi, involves smoothing each segment 

separately and then averaging the smoothed segment estimates. The 

2 

smoothing is accomplished uniformly over each bispectral location’s R 
neighboring frequency points [Ref. 12]: 

K R/2 R/2 

Rl(a)p©2) “ ~2 H 

,-l R r--R/2 s~-R/2 


2nf 2nf, 

where = (“ 7 —(—^)A,,,andR^L. 

JL/ 

In the second approach, denoted ^ 2 ^ unsmoothed segment estimates are 
first coherently averaged: 

K 

(2.19) 


1 


R2(o>i.m2) = 

i - 1 


14 





Then a frequency domain smoothing window is applied through a two- 
dimensional convolution. [Ref. 14] 

Hi-Spec software [Ref. 14] is used in this study to compute the bispectrum by the 
direct method. This program uses the second approach described above, based upon 
(2.19) to reduce the variance of the final bispecirum estimate. There are a few options 
available in this software package to accomplish the frequency domain smoothing. The 
option chosen ihost frequently is the Rao-Gabr window defined by: 


W(m,n) = 


m^ + n^ + mn 


for (m, n) £ G; 


where Af = , or half the Fourier transform length. The set G is the collection of (m, n) 

A 


points satisfying 


2.2. ^ winlen 

m +n +mn^ -^. 


The parameter winlen is the desired window length. Spatial smoothing is applied over 
the (winlen) ^ adjacent frequency points around each bispectrum point. [Refs. 12,14] 


4. Comparison of Methods 

Both the direct method and the indirect method lead to asymptotically unbiased 
estimators, i.e., E {S ^ indirect method has an asymptotic 

variance given by: 

var {Re R (©j, © 2 ) ) = var {Im B (©j, © 2 ) } 

= ^^F(©j)F(®2)F(©i + ©2): 

M K 

where is the energy of the tri-awrelation domain window, P (©) is the true power 
spectrum, Af = 2L + 1 is the Fourier transform size (i.e., m MxM FFT), and is the 
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number of segments averaged [Ref. 1]. The asymptotic variance of the direct method 
when using the smoothing approach described by (2.18) is: 

var {Re (cOp ( 1 ) 2 ) } = var (Im 5 0 ) 2 ) } 

= (o)^)/> (co^)(©1 + ©2); 

“iKLi 

where N is the length of the data sequence, is the number of segments, and L is the size 
of each segment. [Ref. 12] 

The direct method and the indirect method are identical when neither uses a 
window, and the indirect method’s tri-correlation estimate is computed for a number of lags 
equal to the fuU segment length so that L = M - 1 in (2.16). 

C. THE 1-1/2 D INSTANTANEOUS POWER SPECTRUM 

The 1-1/2 D instantaneous power spectrum (1-1/2 Djps) is a combination of the 
standard 1-1/2 D spectrum (1-1/2 Djtd) and the Instantaneous Power Spectrum (IPS) [Refs. 
3,5,16]. It is shown in [Ref. 5] that 1-1/2 Djpj performs better than the conventional 
spectrogram in some respects. When used to process dynamic signals, 1-1/2 D^p^ is 
observed to have a quicker spectral rise time and a quicker fall-off time than the 
spectrogram. The 1-1/2 Djpj method also appears to be good at detecting low SNR 
stationary signals. [Ref. 5] 

The standard 1-1/2 D spectrum is a degenerate form of the bispectrum. When the 
first lag in (2.7) to set to zero, the third order cumulant expression becomes 

C^^^0,/2) = E{A*(n)X(/i)X(n + /2)} (2.20) 

= E{\Xin)\hin + l^)}. 

A one-dimensional Fourier transform of (2.20) produces the 1-1/2 Djjjj spectrum. [Ref. 3] 
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The Instantaneous Power Spectrum is defined as the average of the derivatives of 
Page’s running spectrum, 

,2 

( 2 . 21 ) 




I s (t) exp (-jlnft) dx 
.00 

and Levin’s backward running spectrum. 




(x) exp (-j2nfx) dx 


( 2 . 22 ) 


A discrete version of IPS is obtained by taking the Fourier transform of windowed 
correlation estimates formed using simple lag products vice a sum of lagged products: 

IPS (n, ©) = 2 E (") 2^ (« - ^') + (n) (« + it) } w (it) exp Hmk ); (2.23) 

where w(k) is the window function, and N is the length of the sampled data sequence 
{2f(0),X(l),...,X(iV-l)} [Refs.5,16]. 

The 1-1/2 Dips algorithm is created by substituting the estimate of (2.20) into (2.23) so 
that 

1-1/2 DipJn,<D) = 2 {\X{n)\^)t {n-k) •¥\X{n)fxin + k)}wik)expH(ak) 

(2.24) 


with variables defined as in the IPS equation. [Ref. 5] 


D. INSTANTANEOUS HIGHER-ORDER MOMENT SLICE 

The time-frequency representation introduced in [Ref. 6] is based upon a one¬ 
dimensional Fourier transform of a higher-order moment slice. While satisfying certain 
constraints, a desired number of different moment slice sets are generated in order to form 
an equal number of different but related time-frequency surfaces. As a consequence of 
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the moment slice construction constraints, signal dependent auto-terms are enhanced and 
cross-terms are de-emphasized when the time-frequency surfaces are coherently 
summed. To further control the contrast between auto-terms and cross-terms, each 
moment slice is multiplied by a kernel function before being Fourier transformed. The 
kernel function controls the degree of auto-term and cross-term concentration on the time- 
frequency surface. [Ref. 6] 

The instantaneous higher-order moment slice (MOMS) method is used to detect 
multicomponent linear chirp signals having different instantaneous frequency slopes, or 
slew rates. One characteristic of this time-frequency representation is that auto-term lines 
pass through the origin of the time-frequency plane. The number of auto-terms so 
identified indicate the number of chirps in the signal. In addition, the angle that the 
auto-term line makes with the time axis is related to the instantaneous frequency of the 
chirp. [Ref. 6] 

The IHOMS method is outlined for a received signal, (t), given by 

M 

+ (2.25) 

<■ - 1 

where M is the number of chirps in the received signal. The instantaneous frequency of 
the i*^ chirp is given by: 

/,(')= (2.26) 

= ±[«, + 2V]. 

A general expression for a sampled version of (2.25) is: 

M p 2 ■ 

x(/i) = J^exp{j2n (^)u^+(p }; (2.27) 
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where the sampling period is assumed to be unity, and N is the number of discrete samples 
in the data record. [Ref. 6] 

The following procedure is used to fonn the IHOMS time-frequency surface [Ref. 6]: 


1. The desired order of the moment slice, r, is chosen first. The order should be 
greater than three so that a variety of related time-frequency surfaces can be 
summed in order to form the composite surface. In general, auto-terms stand 
out better when more surfaces are summed. However, if the order is chosen 
too large, cross-term effects increase and may cause spurious peaks in the com¬ 
posite time-frequency surface. 


2. The parameter vector, g 




1 af 


, is generated somewhat arbi¬ 


trarily but subject to the following constraints: 


r-l 


= 1 . 


/-I 


(2.28) 


and 


r-l 


L i""'’ 
/-1 



(2.29) 


3. The moment slice is then calculated for a given instant of time (0, a given 
parameter vector, and a desired number of lags (t): 

= X* (T)x(r + a^T) ..,A:(t+ a^_ jt) . (2.30) 

4. The parameter g in the kernel function, exp , is chosen very small so that 

auto-terms concentrate along straight lines. However, if g is made too small, 
cross-terms show up everywhere on the time-frequency plane and tend to add up 
like auto-terms when the surfaces are later summed to form the composite 
surface. 
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5. TTie Fourier transform of the moment slice is computed and multiplied by the 
kernel function, 

= ^expi-gT^)m^^(x’,t,a^^^)expi-j2nfx)dx. (2.31) 

_oo 

6. Steps three through five are repeated for a different time instant, the same 
parameter vector, the same g, and the same number of lags. Stacking the 
results of (2.31) for successive time instants forms a complex valued time-fre¬ 
quency surface based upon a particular parameter vector. 

7. A new parameter vector is generated and then steps three through six are 
repeated to create another time-frequency surface. This is done again and 
again until the desired number of complex valued surfaces are computed. The 
composite surface is obtained by taking the magnitude of the coherently 
summed parameter vector dependent time-frequency surfaces, 

K 

Git,f) = £Z(r./;fl^*^) . (2.32) 

kml 

where K is the number of different parameter vectors used to generate an equal 
number of related time-frequency surfaces. 

The maxima of the function, 

R2 

Z?(0) = ^ G(vcos0, vsin0), (2,33) 

vtif, 

where 0 is varied between zero and 7i/2 radians, and the radial distance v is varied 
inCTementally between and f? 2 , correspond to chirps. The value of b in the slope term 
of a given chirp’s instantaneous frequency expression, (2.26), can be found from the 0 
producing the associated maximum value in (2.33) using the relation: 
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(2.34) 


bj = jgtanS.; 

where the subscript j denotes a particular chirp, p is the number of lags used to compute 
the moment slice in (2.30), and N is defined as in (2.27). 
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III. PROCESSING GAIN COMPARISONS 


In this chapter we compare the performance of 1-1/2 ^ips to that of the periodogram 
and its time-frequency representation, the spectrogram. The next chapter will examine 
bispectral as well as IHOMS methods, and compare them to 1-1/2 Djpg and the spectrogram 


where appropriate. 

Matlab programs used to conduct the computer simulations for this chapter are 
described first. Theoretical calculations and simulation results are then presented to show 
how each method treats noise and signal components separately. Finally, relative signal- 
to-noise ratio processing gains are estimated. 


A. SIMULATION PROGRAMS 


The 1-1/2 Dips algorithm was originally implemented in [Ref, 5] with the extrinsic 

Matlab function ONE_HALF. Modifications were subsequently made to the program to 
remove the mean from each data segment and to delete some unnecessary smoothing 
steps. The revised ONE_HALF function is contained in Appendix C. The user specifies 
an input data sequence, the size of the data window, a step increment, and a window 
function. ONE_HALF returns the time-frequency representation of the input 
sequence. The ONE_HALF algorithm is based upon (2.24) but is constructed in such a 
way as to avoid imdesirable product term interferences that may form spurious peaks in the 
time-frequency representation. The following sequence is first generated: 


prod\ 


element (0) element (1) 


...element^ 


winlen 

~~T~ 



where. 


(3.1) 
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element(k) = |;e(«)f {jc* (n-^)+ a:(« + ^) } , (3.2) 

winlen is the specified data window length, and jcfn) is the input data segment. A second 
sequence, prod2, is formed by deleting the first element in prodl and then reversing the 
order of the conjugated remaining elements; 


prodl = 


.winlen. .winlen .. 

prodl (— ).pro^l (——^-1). prodl (2) 


(3.3) 


The two sequences are then concatenated and a zero is inserted between them to form a 
proper correlation function having a real valued Fourier transform. The 1-1/2 Dips 
estimate can then be obtained by multiplying the real part of the transformed sequence by 
one-half. However, since the magnitude squared estimate is required later in this study 
to compare the 1-1/2 Dips spectrogram representations, ONE_HAU computes the 
magnitude vice the real part of the transform; 


1 - 1/2 = 


{ \prodl 0 prod '^} 


(3.4) 


where 3 denotes Fourier transform. 

SPECTRO is the extrinsic Matlab function included in Appendix C that implements 
the spectrogram. It has the same input and output format as ONE_HALF. The two 
functions are equivalent with respect to zero padding and transform length so that fair 
comparisons can be made between the two methods. At each step the function computes 
the periodogram; 




winlen —I ^ 

x{n)w{n)exp{-j(Sin) ; 

n-0 


(3.5) 


where the data segment, x(n), is padded with zeros to the transform length of winlen if 
necessary, and w(n) is an appropriate window fimction. A rectangular window (i.e., 
w(n) - 1 over the support of n) is assumed in the estimate variance calculations to follow. 
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B. NOISE ONLY PERFORMANCE COMPARISON 

The variance of the 1-1/2 Djps estimate is derived for a zero mean white Gaussian 
input. The details of this calculation are contained in Appendix A. If the input is real, 
the theoretical variance based on (2.24) is given by, 

var {1-1/2 Dips) =- > (5-6) 

where winlen is the length of the data window and is the variance of the input. For a 
complex input the variance is 

var * [winlen + 2] ia^)^. (3.7) 

The variance of the periodogram as defined by (3.5) with a rectangular window is [Refs. 

2 , 11 ]; 

var {Per} m winlen (a^) . (3.8) 

Both estimate variances are dependent upon the length of the data window and the input 
variance. For a given input variance and increasing window length, the periodogram’s 
variance increases more rapidly than does the estimate variance of 1-1/2 because the 

periodogram’s variance is proportional winler? rather than winlen. However, since the 
1-1/2 Djpj variance is a function of the input variance cubed and the periodogram’s vari¬ 
ance depends upon input variance squared, 1-1/2 Dips’s variance increases quicker than 
the periodogram’s for a fixed window length and increasing input variance. 

Computer simulations confirm the theoretical variance expressions (3.6) and (3.8), 
and the expected trends as window length and input variance are changed. Real valued 
white Gaussian noise is processed by ONE_HALF employing a rectangular window and a 
step increment set equal to the window length so that no overlapping of data segments 
occurs. Since SPECTRO is written for a fixed step increment of one, the periodogram is 
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estimated directly from (3.5) in the Matlab script file used to perform the simulations. 
Simulation parameters, measured estimate variances, and theoretically expected variances 
are shown in Table 3.1. Measured input variances are used to compute the theoretical 
variances. There is good agreement between theoretical and expected values. The 
methods cannot be compared number for number to each other because the periodogram 
variance is a magnitude squared quantity and 1-1/2 Djpj is a magnitude quantity. 
However, comparisons between the methods can be made regarding their relative 
•responses to changing window length and input variance. For a window length of 256, 
1-1/2 Djpj's variance increases by a factor of about 342 when input variance is increased 
from one to seven. For the same variance change, the periodogram’s variance increases 
by a considerably smaller factor of about 49. When input variance is three and window 
length increases from 128 to 512 the periodogram’s variance increases by a factor of 15 
while 1-1/2 Djps only increases by a factor of three. 


TABLE 3.1: THEORETICAL VS. MEASURED 1-1/2 Djps AND PERIODOGRAM 
VARIANCES FOR VARIOUS WINDOW LENGTHS AND INPUT VARIANCES 


input 

variance 

window 

length 

1-1/2 Dips 

theoretical 

variance 

1-1/2 Dips 

measur^ 

variance 

periodogram 

theoretical 

variance 

periodogram 

measured 

variance 

1 

128 

204.3 

205.5 

20,456 

19,487 

1 

256 

394.2 

377.8 

77,549 

93,217 

1 

512 

779.8 

707.6 

303,490 

297,930 

3 

128 

5,515.8 

5,549.5 

184,100 

175,380 

3 

256 

10,644 

10,201 

697,940 

838,960 

3 

512 

21,056 

19,107 

2,731,400 

2,681,300 

7 

128 

70,071 

70,499 

1,002,300 

954,850 

7 

256 

135,220 

129,590 

3,799,900 

4,567,600 

7 

512 

267,490 

242,720 

14,871,000 

14,598,000 
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The 1-1/2 Dipg and spectrogram time-frequency representations of a white Gaussian 
noise input are shown in Figure 3.1. ONE_HALF generates the 1-1/2 D^pj estimate and 
SPECTRO produces the spectrogram. Both representations are formed by applying a 
rectangular window to a single 128 noise sample realization using a 64 point window length 
and a step increment of one. The most striking difference between the two surfaces is that 
1-1/2 Dips consists of several short time duration spikes while the spectrogram exhibits a 
more continuous and smoother surface. The spikier appearance of the 1-1/2 Dips surface 
indicates that the 1-1/2 Dips estimate de-correlates faster than the specirogram. As a 
consequence, the typical 1-1/2 Djps frequency bin contains more independent segments 

than does the typical spectrogram bin. This is significant because a larger reduction in 
variance is realizable from a greater number of independent segments when the power in a 
bin is summed. 

The number of independent segments in a bin is estimated by applying a relation used 
to determine the degrees of freedom for a chi-squared distribution. Given a collection of 
n normally distributed zero mean unit variance random variables, the degrees of freedom 
(DF) represent the number of independent random variables contained in the 
collection. Degrees of freedom are found from the mean and variance of the sum of the 
random variables, s = +X 2 +...using the relation: 

DF = (3.9) 

where and are the mean and variance respectively, of s. [Ref. 7] 

Ten spectrogram and 1 -1/2 Djps representations for different realizations of real white 
unit variance Gaussian noise are used to estimate the degrees of freedom with (3.9). The 
degrees of freedom are calculated for the individual time-frequency cell (i.e., the estimate’s 
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value in one frequency bin at one time step), as well as for the frequency bin. The 
degrees of freedom for a cell is based upon the mean and variance of all the values in the 
frequency bin. The degrees of freedom for a bin is based upon the mean and variance of 
the average frequency bin values for all the bins in the surface. The final degrees of 
freedom estimates are obtained by averaging the degrees of freedom computed for each of 
the ten realizations. Table 3.2 summarizes the results of the simulations for various 
window lengths and window functions. In all cases, SPECTRO and ONE_HALF use a 
step increment of one, and the data sequence is twice the size of the window length. The 
degrees of freedom for a cell in the 1-1/2 Djpj representation is approximately 0.6 while 
the degrees of freedom for the spectrogram cell is about four. The degrees of freedom 
for a bin depends upon the window function. With a rectangular window, a 1-1/2 Djpj 
bin has about 17.5 degrees of freedom and a spectrogram bin has almost six degrees of 
freedom. If a Hamming window is applied to the data, the spectrogram and 1-1/2 Djpg 
bin degrees of freedom figures increase to about eight and 27, respectively. 


TABLE 3.2: DEGREES OF FREEDOM FOR 1-1/2 Dips AND SPECTROGRAM 

ESTIMATES 


winlen 

window 

function 

ONEHALF 

DFceU 

ONEHALF 

^F^in 

SPECTRO 

DFceU 

SPECTRO 

DFbin 

64 

rectangular 

0.6582 

16.3737 

4.1109 

5.8239 

64 

Hamming 

0.6185 

26.9717 

4.4057 

7.7394 

128 

rectangular 

0.6326 

18.7693 

4.0822 

5.9613 

128 

Hamming 

0.5974 

27.5028 

4.1742 

8.0081 


The data sequence length divided by the bin degrees of freedom indicate the size of 
the independent segments in a typical frequency bin. For example, the size of an 
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independent segment in the 1-1/2 Djpg estimate based on a window length of 64, a 128 point 

sequence, and a rectangular window is 128/16.3737 = 8, For the same parameters, the 
spectrogram independent bin segment size is about 22. These figures seem reasonable 
based upon “eyeball” estimates of independent bin segment size in Figure 3.1 contour plots. 

When an estimate’s frequency bin values are averaged, the anticipated variance 
reduction factor is given by the ratio of that estimate’s DF^in to DF^^u values. The 
expected variance reduction factor for 1-1/2 Dips about 29 for a rectangular window and 
45 for a Hamming window. The spectrogram has an expected variance reduction factor 
of 1.5 for a rectangular window and two for a Hamming window. 

Computer simulations using ONE_HALF and SPECTRO confirm that an exploitable 
difference in variance reduction capability exists between 1-1/2 Dip^ and the spectrogram 
when frequency bin averages of a time-frequency surface are computed. Gaussian white 
noise is simulated to serve as a noise only, real input signal. Both methods employ a 
rectangular window on the same input signal. The step increment is specified to be one 
for each method and the data window length is consistently set equal to half of the data 
sequence length. For the first set of simulations, noise variance is fixed while the 
transform length is varied. The simulations require that each method generate a time- 
frequency representation for fifteen different input sequences. The fifteen representations 
are then averaged to form one time-frequency surface. For this, and all remaining 
simulations performed in this study, the averaged 1-1/2 Djps surface is squared so that its 
amplitude can be fairly compared to the spectrogram’s amplitude. As a consequence of 
computer program design, data sequence length, and choice of step size the first and last 
winlenll time steps are calculated using zero-padded data segments. These time steps are 
deleted and the average value of each frequency bin on the average surface is calculated. 
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Figures 3.2 through 3.6 show these average frequency bin values for an input noise 
variance of one and window lengths of 32, 64, 128,256, and 512, respectively. The bin 
corresponding to dc. and the bin corresponding to the fold-over frequency are omitted from 
the plots. To allow for a meaningful comparison, the mean frequency bin average is 
subtracted from each plot. Examination of Figures 3.2 through 3.6 indicate that 1-1/2 Djps 
exhibits a lesser output variation than does the spectrogram as the transform length 
increases. TWs observation is confirmed in Figure 3.7 which plots the measured output 
variance versus transform length for the first simulation set. 


^-^/Z Dips 



Spectrogram 



Figure 3.2: Averaged output variance for winlea - 32 and = 1 • 
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Figure 3.7: Spectrogram (dashed) and 1-1/2 Djpg (solid), output variance vs. window length 

foraf^ = l. 


The simulation was repeated for different input noise variances and a fixed window 
length. Figure 3.8 shows the relationship between output variance and input variance for 
a window length of 128. The point at which the two curves in Figure 3.8 intersect can be 
considered a performance breakpoint. At input variances below the breakpoint value (i.e., 
about 1.5 from Figure 3.8), the 1-1/2 estimate has less variation than the 
spectrogram. When the input variance is larger than the breakpoint value the spectrogram 
exhibits less variation than 1-1/2 D^p^, Other simulations show that breakpoints for 
window lengths ranging from 32 to 512, although gradually increasing as window length 
increases, all correspond to an input noise variance of approximately 1.5. 

Figure 3.8 shows that the spectrogram has a variance of about 1000 when the input 
variance is about 1.3. Considering the processing scheme, the expected variance is the 
theoretical estimate variance given by (3.8) reduced by two factors. The first factor is 
equal to the number of time-frequency surfaces averaged to obtain a representative surface. 
The second factor is the degrees of freedom based variance reduction factor estimated 
earlier. Therefore, the expected spectrogram output variance is given by: 
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2 'I 2 

winlen 

var ISpectroym-—^^ 

_ 1282- 1.32 
15-1.5 
= 1231; 

where the numerator is (3.8), nsurf is the number of time-frequency surfaces averaged in 
the simulation, and VRF is the degrees of freedom variance reduction factor for a 
rectangular window. The expected and actual variances for the spectrogram agree fairly 
well considering that the comparison is based upon approximate values. 

The same comparison is made for 1-1/2 Dip^. According to Figure 3.8, an input 
variance of approximately 1.45 results in a 1-1/2 Djpj estimate variance of about 
1000. The theoretical variance given by (3.6) is reduced by the same factor of 15 used in 
the spectrogram calculation since this factor depends only upon the number of time- 
frequency surfaces averaged. The theoretical variance is also reduced by a degrees of 
freedom variance reduction factor that is larger than the spectrogram’s (29 vice 1.5). To 
account for squaring the output of ONE_HALF in the simulations, (3.6) is also 
squared. The expected 1-1/2 Djp^ variance based on the processing scheme is: 


var {l-mOips} s 


[54 + Swinlen] ^ (a^) ^ 

4^ • nsurf - VRF 
[54-h 6 - 128 ] 2 ( 1 . 45 )® 
4^•15 ■ 29 


= 902. 


The expected variance agrees fairly well with the simulation variance for the 1-1/2 Djpj 
estimate considering that approximate values are being compared. 
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Figure 3.8; Spectrogram (dashed) and 1-1/2 (solid), output variance vs. input variance 

for winlen = 128. 


C. SIGNAL ONLY PERFORMANCE COMPARISON 

Appendix B details the calculations made using (2.24) to determine the theoretical 
1-1/2 Djpj estimate for a signal only. For a real valued sinusoid given by 

x{n) - cos(27t^«), (3.10) 

the theoretical estimate is: 

N^tn 

1-1/2 Djpj,(n,a>) = 2^^®® (27i-n) {6(Q)-m)+6[£0-(N-m)]}. (3.11) 

The periodogram is defined in (3.5). For the same signal, the theoretical value of the 
periodogram is the magnitude squared of the Fourier transform of the input signal: 

N 

Per(n,ai) = (^) {5(<o-m)+6[co-(N-m)] } ( 3 . 12 ) 

If the signal consists of a single complex sinusoid, 

x(n) =exp(J2n~n), (3.13) 
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the theoretical 1-1/2 Djpj estimate is given by: 

1-1/2 Djpj, (n,<») = iVcos (271^/1)6 (©-m); (3.14) 

and the theoretical value of the periodogram is: 

Per(n,©) s iV^6(©-m) (3.15) 

Results of computer simulations using SPECTRO and ONE_HALF validate the 
theoretical calculations for the periodogram and 1-1/2 Dipg, respectively. For the 

simulations, a 128 point data window is stepped through a 256 point data sequence one 
point at a time. Both methods apply a rectangular window function to the data. The first 
and last winlen/2 time steps are discarded before computing any averages or plotting any 
results. The 1-1/2 Djps estimate is squared to allow a fair comparison with the 
spectrogram. Figure 3.9 shows the mesh plot and corresponding contour plot of the 
1-1/2 Djpj estimate for the real signal. The sinusoidal signal is centered in frequency bin 
20. Bin 20 and the average of each frequency bin are plotted in Figure 3.10(a) and 
Figure 3.10(b), respectively. The maximum power at any one time step is approximately 
4042. This agrees fairly well with the magnitude squared value of (3.11) which is equal 
to half of the window length squared, or 4096. The intra-ridge modulation exhibited by 
bin 20 in Figures 3.9 and 3.10(a) results from the cosine cubed modulating term in (3.11). 
The average power of bin 20 is approximately 1280 as shown in Figure 3.10(b), 

Figures 3.11 and 3.12 use the same arrangement of plots to display spectrogram results 
for the real signal. The maximum power agrees exactly with the theoretical value of4096 
calculated using (3.12). Figure 3.12 shows that the average power of bin 20 is equal to 
the maximum power. The maximum power of each method is about the same, however, 
the average power of 1-1/2 Djpj is approximately one-third as much as the spectrogram’s 
average power due to the intra-ridge modulation effect. 
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(a) (b) 

Figure 3.9:1-1/2 Djpj time-frequency representation of a real sinusoid. 




(b) 


Figure 3.10: 1-1/2 Djpj signal power in bin 20 (real sinusoid). 




































































































































































































When the same simulations are run with a complex sinusoidal input, ONE_HALF 
produces a maximum signal power approximately equal to the length of the data window 
squared which agrees with the squared magnitude of (3.14), As a result, ONE_HALF’s 
maximum power is about equal to the spectrogram’s for both real valued and complex 
valued signals. However, for the complex valued sinusoid die average power in bin 20 of 
the 1-1/2 Djpj estimate is about half that of the spectrogram’s rather than a third as it was 
for the real sinusoid. This is a consequence of the different modulating terms in (3.11) and 
(3.14). The real signal is modulated by a cosine cubed function which has a magnitude 
squared average value of about 0.345 for a unit amplitude cosine over one period. In 
contrast, the complex valued signal is modulated by a cosine function having a magnitude 
squared average value of about 0.52, 


D. SNR PROCESSING GAIN COMPARISONS 


The ability of a technique to detect a signal in noise is measured in two ways. The 
first measure is a processor signal-to-noise ratio (PSNR), or simply the ratio of processor 
signal power to processor noise power: 


PSNR = 



(3.16) 


where Pn is the noise variance. The second approach, denoted PDIF, measures the power 
of the signal in terms of the number of noise standard deviations above the mean noise 
power level: 


PDIF = —~ 


(3.17) 


where p, is the processed noise mean power and is the standard deviation of the 


processor output noise power. The ratio of either (3.16) or (3,17) for two processing 
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techniques forms a figure of merit used to compare methods. To compare 1-1/2 Djpjj and 
the spectrogram using (3.16) the processing gain ratio is defined as; 


_ ^^^^ONEHALF 

PGsm - psnRspectro 


(3.18) 


When (3.17) is the measure of interest, the processing gain ratio is given by; 


P^^^ONEHALF 


"DIF pQif 


(3.19) 


SPECTRO 


Computer simulations are used to measure the processing gain for various window 
lengths and input noise variances. The input consists of a fixed amplitude real sinusoid 
centered in frequency bin 20 and mixed with Gaussian white noise. For this scenario the 
input signal-to-noise ratio (SNR) is given by; 

SNR = I01og|^^j (3.20) 


Fifteen realizations of each simulation are averaged in order to obtain representative 
results. As before, the input data sequence is twice the length of the data window, a 
rectangular window is applied to the data, and the data window is stepped through the 
input sequence one point at a time. Figure 3.13 shows the resulting plots for a data 
window length of 128, and an input SNR of -3dB. The fifteen-realization averaged 
1-1/2 Dips time-frequency surface is shown in Figure 3.13(a) and its corresponding 
frequency bin average is displayed in Figure 3.13(b). Figures 3.13(c) and 3.13(d) are the 
averaged spectrogram time-frequency surface, and its frequency bin averages, 
respectively. Figures 3.14, 3.15, 3.16, 3.17, and 3.18 display the results for the same 
window length but for different input SNRs of about -lOdB, -15dB, -17dB, -18.5 dB, and 
-19.5 dB, respectively. 
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Processing gain comparisons using (3,18) and (3.19) were made of the frequency bin 
averages as a function of SNR for the signals shown in Figures 3.13 through 3.18. The 
average value of the signal bm was used for P^ig to compute PSNR and PDIF. The mean, 
standard deviation, and variance of the noise were computed using all bins except the one 
containing the signal. Table 3.3 summarizes the results of these processing gain 
computations. 

The PGsnr performance measure is consistent with results obtained from the noise 
only simulations and the signal only simulations. As expected from these earlier 
simulations, POg^R is greater than one (indicating that 1-1/2 ^ips performs better than the 
spectrogram) when the input variance is smaller than about 1.5. This is agrees with the 
breakpoint determined from the analysis of Figure 3.8 earlier. Based upon the time- 
frequency representations however, PGuip appears to be a more useful indicator of 
detection performance over a much wider range of input SNRs. 

The PGdip column of Table 3.3 shows that the 1-1/2 Djps frequency bin average signal 
peak estimate, when measured in terms of noise standard deviations above the mean noise 
level, is about 1.6 times higher than it is for the spectrogram estimate when the input SNR 
is about-3dB. As input SNR decreases the relative processing gain also decreases. Both 
methods are equivalent when input SNR is about -19dB. The spectrogram becomes the 
better estimate by this standard at SNRs lower than -19dB for this window length. 
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(C) (d) 

Figure 3.16; Sinusoid in bin 20 with SNR s -17dB; (a) 1-1/2 Djp^ time-frequency 
representation, (b) 1-1/2 Djpj bin averages, (c) spectrogram time-frequency representation, 

and (d) spectrogram bin averages. 
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(C) (d) 

Figure 3.17: Sinusoid in bin 20 with SNR s -18.5dB; (a) 1-1/2 Djpj time-frequency 
representation, (b) 1-1/2 Djpj bin averages, (c) spectrogram time-frequency representation, 

and (d) spectrogram bin averages. 



























X 10 































TABLE 3.3: RELATIVE PROCESSING GAINS FOR VARYING INPUT SNR AND 
A FIXED DATA WINDOW LENGTH OF 128 


input SNR 
(-dB) 

PGsnr 

POIFonehalf 

PDIFspectro 

PGdif 

corresponding 
figure no. 

3.01 

1.4270 

217.8826 

132.4518 

1.6450 

3.12 

10.00 

< 1 

37.9275 

29.3067 

1.2942 

3.13 

14.77 

< 1 

12.8635 

11.2700 

1.1414 

3.14 

16.99 

< 1 

8.0222 

7.4846 

1.0718 

3.15 

18.45 

<1 

5.9961 

5.8097 

... 

1.0269 

3.16 

19.54 

<1 

4.8344 

4.8552 

0.9957 

3.17 


Results of signal-in-noise simulations run for a smaller 64 point data window and a 
larger 256 point data window are consistent with those obtained for the 128 point 
window. However, the range of input SNRs over which the 1-1/2 Dips method yields a 
larger processing gain than does the spectrogram depends upon window length. For the 
64 point window, 1-1/2 Djpj yields a larger processing gain than the spectrogram only for 
input SNRs greater than about -12.5dB. As already demonstrated, for a 128 point window 
1-1/2 Dips ^ larger processing gain than the spectrogram for input SNRs greater than 
about -19dB. Simulations using a 256 point window showed that 1-1/2 Dips performs 
better than the spectrogram in terms of processing gain even when the input SNR is less 
than -21dB and the signal cannot be readily discerned from the noise in the averaged 
frequency bin plots. In effect, for this type of signal and processing scheme using a 256 
point data window, 1-1/2 Djps's processing gain is between 1.15 and 1.64 times larger than 
the spectrogram’s for input SNRs between -2 IdB and -3dB. The simulations also indicate 
that the processing gain performance of 1-1/2 Djpj improves relative to the spectrogram 
with increasing window length. 
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IV. SIMULATION RESULTS AND ANALYSIS 


Simulation results are presented to show how 1-1/2 D^pg, the bispectrum methods, 
IHOMS, and the spectrogram deal with three types of signals. The spectrogram, 
bispectrum, and 1-1/2 D^pg methods are first compared using a signal consisting of multiple 
unrelated stationary sinusoids that are mixed with white Gaussian noise. Next, the 
bispectral representation of noise free signals containing harmonically related stationary 
components is studied. Lastly, a received signal containing multiple linear FM 
components mixed with white Gaussian noise is processed using IHOMS, 1-1/2 Djpg, and 
the spectrogram. , 

A. STATIONARY SINUSOIDS 

The spectrogram is computed using the Matlab function SPECTRO, and the 1-1/2 D^pg 
method is implemented via ONE_HALF. Their time-frequency representations for a 
single sinusoidal signal are discussed in Chapter IE. 

The bispectrum methods produce a frequency-frequency representation instead of a 
time-frequency representation. The direct method for computing the bispectrum is 
realized through the Hi-Spec Matlab function BISPEC„D of [Ref. 14]. This function 
returns the complex bispectrum matrix and a frequency axis labeling vector. The user 
specifies an input data sequence, the desired Fast Fourier Transform (FFT) size, the number 
of samples per segment, the amount of overlap between data segments, and the desired 
form of the frequency smoothing window. The extrinsic Matlab function INDBIS 
contained in Appendix C implements the indirect method of calculating the bispectrum. It 
returns the complex bispectrum of the signal supplied by the user. The FFT size, the 
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number of data samples per segment, and the desired tri-correlation window function are 
specified by the user. The three window options available are the unit hexagonal window, 
the Parzen window, and the optimum window (also known as the minimum bias supremum 
window). 

Both bispectrum methods are used to generate bispectral representations of a noise 
free sinusoid located at frequency bin 20. The input sequence consists of 256 
samples. Both methods apply a 128 point FFT to each of four non-overlapping data 
segments containing 64 data points each. The indirect method representations are formed 
from tri-correlation functions based upon 31 lags. Figure 4.1 shows the full bispectrum 
plane representation of a real sinusoid generated by the direct method without a 
window The sinusoid is indicated by the peak in the first quadrant located at 
(kj, ^ 2 ^ = (20,20), and associated symmetry peaks like the one in the third quadrant at 

(ki,/:2) = (-20,-20). 

The representation of a real sinusoid at bin 20, formed by tiie indirect method using a 
unit hexagonal window is displayed in Figure 4.2. The rectangular window used by both 
ONE_HALF and SPECTRO, the unit hexagonal window used by INDBIS, and the 
unwindowed BISPEC_D are essentially equivalent window functions so that 
representations formed using them are considered unwindowed. Two closely spaced 
peaks are discernible at the signal locations in the indirect method’s representation. This 
is a characteristic of this method which depends upon the number of lags used to compute 
the tri-coirelation function. The two peaks at each location become less distinct as the 
number of lags used to compute the tri-correlation sequence is increased. Another 
difference between Figure 4.1 and 4.2 is that the indirect representation exhibits peaks at 
(kj, ^ 2 ) = (20, -20) and (kj^, k^) = (-20,20) which are not noticeable in the direct 
method’s representation. It turns out that both of the methods are sensitive to the phase 
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of a real signal, and that the methods produce slightly different representations for the same 
phase. The signal used to create the bispectrum in Figure 4.1 has a different phase than 
the signal which the bispectrum in Figure 4.2 is based upon. Both phase values are 
specially chosen to chosen to minimize any peaks on the = 0 and k 2 = 0 axes. In the 
case of the direct method’s representation, this also minimized the peaks at 

^: 2 ) = (20, -20) and (“20, 20). Bispectral peaks located on either 

the or the ^2 in subsequent plots are referred to as zero axis peaks. 

The unwindowed full bispectrum representations for a complex sinusoid obtained by 
the direct method and indirect method are shown in Figure 4.3 and 4.4, respectively. Like 
the real signal representation, the complex signal produces a peak at (kj, ^ 2 ) = (20,20). 

The complex signal’s bispectrum does not exhibit as many symmetry peaks as does the real 
signal’s representation because the bispectrum of a complex signal is symmetric about only 
one symmetry line and in only one quadrant. Which line, and which quadrant are 
determined by the conjugation scheme used to compute the bispectrum [Ref. 9). For the 
equivalent schemes used by INDBIS and BISPEC_D which are described in Chapter H, 
symmetry exists about the k^ - ^2 line in the first quadrant. The peaks at 
{k^,k-^ = (0,20) and (^^,^ 2 ) “ (20,0) in both representations are zero axis 

peaks. Unlike real signal representations, zero axis peaks in the bispectrum of a complex 
signal have a consistent amplitude regardless of signal phase. 

The first bispectral quadrant is essentially the same for either a real or a complex 
sinusoid, and is sufficient for measuring a sinusoidal signal’s frequency. The only 
additional information not contained in the first quadrant is whether or not the signal is 
complex or real. The remaining bispectrum figures in this thesis will display just the first 
quadrant. 
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figure 4.1: Full bispectrum of a real sinusoid at bin 20 (imwindowed direct method). 


X 10 



eo 


40 


20 


g 
^ o 

<D 
§■ 

-20 


•4-0 


50 


--1 

« 




. % . 

A . 






4 . 

. % 

.«>. 

..... 





:_i _^— 


k2 axis 


-50 O 50 

k1 frequency bln axis 


(a) (b) 

Figure 4.2: Full bispectrum of a real sinusoid at bin 20 (unwindowed indirect method). 
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The next four examples of single real valued sinusoid bispectrums are intended to 
show what windowing can do for each method. At the same time, the effects of signal 
phase on the bispectral representations are noted. Figure 4.5 shows the bispectrum as 

computed by the direct method and smoothed by a Rao-Gabr window over 5 adjacent 
frequency points. The signal is a real sinusoid located at bin 20 with a phase angle of 0.484 
radians. The same processing scheme yields the representation shown in Figure 4.6 if the 
signal phase is changed to 4.3151 radians. The indirect method, employing a Parzen 
window on a real signal with a phase equal to 0.484 radians produces the representation 
shown in Figure 4.7. When the signal with a phase of 4.3151 radians is processed by the 
indirect method using an optimum window (the optimum window and Parzen window are 
not significantly different for this type of signal), it produces the bispectrum shown in 
Figure 4.8. Windowing helped to make the two distinct peaks at * (20,20) 

appear as one in the indirect representation. The zero axis peaks of the indirect 
representation incorrectly appear to be offset from the zero axes because the first quadrant 
plots do not show the adjacent peaks in the second and fourth bispectral quadrants. The 
signal frequency is more accurately and more clearly depicted by the direct method’s peak 
locations. In addition, the direct method is faster than the indirect method as implemented 
in this study. However, if the signal is known to be either complex or real, the indirect 
method can be implemented in such a way as to take advantage of the bispectrum symmetry 
in order to reduce computational costs. This is done in [Ref. 14] to implement the indirect 
method for real signals in Hi-Spec’s Madab function BISPECJ. 

The bispectrum, spectrogram, and 1-1/2 Djps methods are now used to process a signal 
consisting of three complex valued sinusoids in additive white Gaussian noise. The three 
frequencies of the signal are at bins 10.65,28.05, and 54.7. Noise variance is fixed while 
signal amplitude is varied to achieve different input SNRs. The phase angle of each signal 
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Figure 4.5: Rao-Gabr smoothed direct bispectrum of real sinusoid (0.484 radians). 
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Figure 4.6: Rao-Gabr smoothed direct bispectrum of real sinusoid (4.3151 radians) 





































































component is distributed uniformly on the interval (0, 27t). The data sequence has 256 
samples. The 1-1/2 spectrum is computed by ONE_HALF, the spectrogram is 
generated by SPECTRO, and the bispectrum is calculated by the direct method using 
BISPEC_D. A step increment of one is employed for both SPECTRO and ONE_HALE, 
while a 50% overlap is specified for BISPEC_D. To allow for a fair comparison, all 
three methods are unwindowed. Fifteen representations for each method are computed 
using different noise realizations. These surfaces are then averaged to obtain either a 
representative time-frequency surface, or a representative frequency-frequency surface, as 
appropriate. The averaged surface, along with its frequency bin average is displayed for 
various input SNRs and window length. The bispectral frequency bin average is the 
average power in each ki frequency bin. A horizontal line on the frequency bin plot is 
placed three noise standard deviations above the noise mean to provide a reference 
point. The noise statistics are computed after removing outliers from the frequency bin 
average values. Figure 4.9 shows the spectrogram representation for a 128 point window 
length, and an input SNR of -6dB. The 1-1/2 Djpj spectrum, and the bispectrum for the 

same input SNR and window length are displayed in Figures 4.10 and 4.11, respectively. 
The signal components are clearly detected by aU three methods but the symmetry of the 
bispectrum causes extra peaks that extend above the three standard deviation reference 
line in Figure 4.11(c). 

Figures 4.12, 4.13, and 4.14 show the spectrogram, the 1-1/2 Dipj, spectrum and the 
bispectrum, respectively, for an input SNR of-18dB and a window length of 128. At this 
lower SNR it is difficult to distinguish the signal components in the methods’ mesh plots. 
Figure 4.12(a), 4.13(a), and 4.14(a). However, the sinusoids are discernible in the contour 
plots and bin average plots of all three methods without any spurious peaks exceeding the 
three standard deviation reference line. The slight processing gain advantage of the 
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Figure 4.11: Bispectrum of three complex sinusoids at an input SNR of -6dB. The 
representation is formed using a 128 point data window. 









1-1/2 Dips niethod over the spectrogram at this SNR and window length is barely 
noticeable in Figures 4,12(c) and 4.13(c). But, the signal peaks and the larger noise 
peaks in the spectrogram plot are closer to the three standard deviation reference line than 
they are in the 1-1/2 Dips plo*. Relative to the spectrogram, the bispectrum signal peaks 
have the same height above the reference line but the noise peaks are closer to the line. 

Figure 4.15 displays just the frequency bin plots of the three methods for an input SNR 
of -21dB. The spectrogram has a slight advantage over 1-1/2 Dips this SNR in terms of 
the proximity of signal peaks and noise peaks to the three standard deviation line. This is 
consistent with Chapter III simulation results which indicate that the relative processing 
gain favors the spectrogram below -19.5dB for this window length. The bispectrum has 
relatively weaker signal peaks and larger noise peaks than either of the other two methods. 

Figure 4.16 shows the frequency bin average plots for the same SNR of -21dB but a 
larger data window of 256 points. The 1-1/2 Djp^ spectrum, as expected for this larger 
window length, appears to perform slightly better than the spectrogram. However, one 
noise peak at bin 60 in the 1-1/2 Djp^ plot just reaches the three standard deviation line in 
Figure 4.16(b). The bispectrum has strong signal peaks and relatively low noise peaks 
except for the one at about bin 118 that exceeds the reference level. 

Figure 4.17 displays the frequency bin plots for an SNR of -26dB and a window length 
of 256. For all three methods, only the first signal component exceeds the reference level. 
The 1-1/2 Dips representation’s signal peak has a slightly larger relative height above the 
reference line than either of the other methods. No one method distinguishes itself with 
respect to suppressing the noise peaks in this case. 
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The representation is formed using a 128 point data window. 
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(c) 

Figure 4.14: Bispectrum of three complex sinusoids at an input SNR of -18dB. The 
representation is formed using a 128 point data window. 
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Figure 4.15: Frequency bin averages for an input SNR of-2ldB and a 128 point window 
length; (a) spectrogram, (b) 1-1/2 Dip^, and (c) bispectrum. 
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Figure 4.16; Frequency bin averages for an input SNR of -21dB and a 256 point window 
length; (a) spectrogram, (b) 1-1/2 Dipj, and (c) bispectrum. 
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Figure 4,17: Frequency bin averages for an input SNR of -26dB and a 256 point window 
length; (a) spectrogram, (b) 1-1/2 Djps, and (c) bispectrum. 
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B. HARMONICALLY RELATED SINUSOIDS 

One of the main advantages of higher than second order spectra lies in their ability to 
preserve true phase information [Ref. 1]. A common application found in the literature to 
test this ability is the detection of quadratic phase coupling [Refs. 1,3,4,11]. Three 
sinusoids are quadraticaUy phase coupled when they are both harmonically related and 
phase related. Harmonically related means that one of the sinusoids has a frequency equal 
to the stun of the other two. Similarly, phase related means that one of the phases is equal 
to the sum of the other two. [Ref. 1] 

Simulations involving BISPEC_D show that the bispectrum is sensitive to signals 
containing harmonic components in general and not just when quadratic phase coupling 
occurs. A signal consisting of two sinusoids is considered initially. The first sinusoid 
has a frequency corresponding to bin 15, the second sinusoid’s frequency corresponds to 
bin 40. Both sinusoids have uniform random phases independently distributed over the 

interval (0, 2?:). The bispectrum representation is generated using BISPEC_D with a 3^ 
Rao-Gabr window, a 64 point data window, and 128 point Ph is. Figure 4.18 shows the 
bispectrum representation of this signal. Both components of the signal are indicated by 
the vertically aligned pair of peaks in ki bins 15 and 40, and/or by the horizontally aligned 
peaks in ^2 t>ms 15 and 40. Figure 4.19 displays the results for the same simulation when 
the second sinusoid’s frequency is changed so that it now corresponds to bin 30 and 
becomes a harmonic of the first component. Only the one peak at = k 2 ~ 15 actually 
stands out on the bispectrum surface. The bispectrum is blind to the sinusoid at bin 30 
because it is the highest harmonic component of the first sinusoid contained in the signal. 
A signal comprised of three harmonic components located at bins 15, 30, and 45 is now 
processed under two phase related conditions. For the first condition, all the phases are 
i.i.d. uniformly over (0,27t). The bispectrum representation for this signal is shown in 


68 



K2 bln axis 





































Figure 4.20. TTie highest harmonic component is again not represented by a significant 
peak located in bin 45. However, the bispectrum does indicate by its appearance that the 
signal contains harmonic components based on the fundamental frequency corresponding 
to bin 15. If the components were not harmonics of this fundamental fi'equency the 
bispectrum surface would have an appearance similar to that previously shown in Figure 


4.11. 



(a) (b) 


Figure 4.20: Three harmonic sinusoids in bins 15,30, and 45 with imrelated phases. 

Figure 4.21 displays the bispectrum for the three harmonic component signal under the 
second phase condition. Here the third sinusoid has a phase equal to the sum of the other 
two component phases. This situation satisfies the definition of quadratic phase 
coupling. There is essentially no difference between this bispectrum and the one 
generated for the unrelated phase condition. When the signal components are harmonics 
of the same fundamental frequency, the bispectrum representation does not discern 
quadratic phase coupling. 
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(a) (b) 

Figure 4.21: Three quadratically phase coupled harmonic sinusoids in bins 15,30, and 45. 


If the signal components are harmonically related but are not harmonics of the same 
fundamental frequency, the bispectrum does distinguish quadratic phase coupling. This 
frequently considered situation in the bispectrum literature is depicted in Figure 4.22. The 
signal components are quadratically coupled but the signal frequencies are 15, 25, and 40 
rather than harmonics (i.e., 15,30, and 45). This representation has just two peaks located 
at the intersection of the ki and ^2 bins where the first two sinusoids reside. 



(a) (b) 

Figure 4.22; Quadratically phase coupled sinusoids located in bins 15,25, and 45. 


71 

























C. MULTI-COMPONENT LINEAR CfflRPS 


A computer simulation based upon Example 1 of [Ref. 6] is used to compare the 
processing ability of IHOMS, 1-1/2 D^p^, and the spectrogram when the received signal 
consists of two linear FM signals mixed with white Gaussian noise. In terms of the 
variables in (2.27), 


M 

x{n) ^ exp {j2n 
/-I 






M is two, = 8, = 94, 1^2 = 1^6, ^2 ” 30, and N = 256. With a record length 

equal to one second and a sampling rate of 256 samples/second, chirp 1 exhibits a slew 
rate of 188 Hz/second, starts at a frequency of 8/256, and finishes at 196/256. Likewise, 
chirp 2 has a slew rate of 60 Hz/second, starts at a frequency of 196/256, and finishes at a 
frequency equal to 256/256. The parameter n,- is the fractional starting frequency of the 
chirp times N. The parameter &,■ is equal to half of the fractional frequency difference 
of the chirp times N, or simply half of the i* chirp’s slew rate. 

The IHOMS method was applied first by following the procedure outlined on pages 
19 through 21. Choosing to work with a fourth order moment slice of the signal fixes the 
length of the parameter vector, a, to four (i.e., g = a 2 , ). GEN_A is an 

extrinsic Matlab function in Appendix C used to generate a desired number of different 
parameter vectors having elements that satisfy (2.28) and (2.29). Three parameter vectors 
are formed for the first simulation. The remaining steps in the procedure are carried out 
by the extrinsic Matlab function ATH_IMS in Appendix C. One-hundred lags are used to 
compute the fourth-order moment slice by (2.30) and one time-frequency representation is 
calculated in accordance with (2.31) for each of the three parameter vectors. The value 
for the parameter g in the exponent of the kernel function in (2,31) is set equal to 0.0008 
because the same value was found to give acceptable results in Example 1 of [Ref. 6]. The 
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final time-frequency representation is obtained by computing the magnitude of the three 
parameter vector specific surfaces. Figure 4.23(a) shows an IHOMS time-frequency 
representation based upon a fourth order moment slice and three parameter vectors for a 
noise-free signal consisting of chirp 1 and chirp 2. The two chirps are indicated by the 
two contour lines radiating from the origin. The remaining tines on the plot are cross¬ 
terms. Figure 4.23(b) is a plot of the one-dimensional function Z)(0) given by 
(2.33). This function represents the magnitude of the time-frequency surface along radial 
tines passing through the origin of the surface. The first maximum in Figure 4,23(b) 
occurs at 0.09 radians and corresponds to the b parameter of chirp 2 through (2.34): 

2 

Applying (2.34) to the location of the second maximum, 0.28 radians, an estimate of the b 


parameter of the first chirp, hi = 94.23, is obtained. 




(a) 


(b) 


Figure 4.23: Three parameter vector IHOMS representation for a noise-free signal. 
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Figure 4.24 shows the three parameter vector IHOMS representation for the same two 
chirps mixed with white Gaussian noise at an input SNR of -3dB. The chirps cannot be 
distinguished from the background noise. 
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(b) 


Figure 4.24: Three parameter vector MOMS representation when SNR = -3dB. 


Signal detection is improved by summing more parameter vector specific surfaces in 
order to form the MOMS representation. This is accomplished by creating more 
parameter vectors. The MOMS representation in Figure 4.25 uses 21 parameter vectors 
on the two chirp signal at an SNR of -3dB. The chirps are once again discernible. Even 
with 21 parameter vectors the detection capability degrades significantly as the SNR 
becomes worse. Whether the chirps can still be distinguished in Figure 4.26 when the 
SNR is -5.5dB is questionable. At lower SNRs they definitely become lost in the noise. 

The l-l/2Dip5 and spectrogram estimates of the double chirp signal are obtained using 
ONE_HALF and SPECTRO respectively. A rectangular window function, a 32 point 
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'igure 4.25: Twenty-one parameter vector IHOMS representation when SNR = -3dB 
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Figure 4.26: Twenty-one parameter vector IHOMS representation when SNR = -5.5dB 










data window, and a step increment of one are the specified function options used in the 
simulations. Figures 4.27, 4.28, and 4.29 show the spectrogram time-frequency 
representations for the noise-free signal, an SNR of -3dB, and an SNR of -5.5dB, 
respectively. Figures 4.30, 4.31, and 4.32 display the 1-1/2 D^pj representations for the 
noise-free case, an SNR of -3dB, and an SNR of -5.5dB, respectively. Unlike the 
IHOMS representation, both 1-1/2 Djpj and the spectrogram indicate fractional start and 
finish frequencies for each chirp in addition to slew rate. For a given noise level, the 
detection capabilities of 1-1/2 D^pj, the spectrogram, and the 21 parameter vector IHOMS 
are comparable. However, IHOMS has a much higher computational cost than either of 
the other methods. Generating a 21 parameter vector IHOMS surface based upon a 
moment slice calculated from 100 lags requires 2,100 100-point Fourier transforms. By 
comparison, the spectrogram and 1-1/2 Djps require only 256 32-point Fourier transforms 
to generate an equivalent time-frequency representation. 




(a) (b) 


Figure 4.27: Spectrogram of noise-free two chirp signal. 



frequency bln 


frequency bin 



Figure 4.29: Spectrogram of two chirp signal with SNR = -5.5dB 


































Figure 4.32: 1-1/2 Djpg representation of two chiip signal with SNR = -5.5dB. 


IHOMS representations indicate only the slew rate of chiips and not their start/stop 
frequencies. As a consequence, IHOMS cannot distinguish between two or more 
different chirps having the same slew rate. Figure 4.33 displays the IHOMS 
representation for two chirps. The first chirp starts at a frequency of 64/256 and finishes 
at 192/256. The second chirp starts at 128/256 and finishes at 256/256. They both have 
the same slew rate of 128Hz/second. The IHOMS representation detects the presence of 
only one chirp in the signal. The maximum of £>(0) occurs at 0.19 radians which 
corresponds to the proper b parameter of 64. The spectrogram and 1-1/2 Djpg 
representations indicate the presence of both chirps. The 1-1/2 Djpg time-frequency 
surface of this signal is displayed in Figure 4.34. 
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The IHOMS algorithm also requires a more intricate sampling scheme than does the 
spectrogram and 1-1/2 Djpj methods. The received signal has to be sampled at a much 
higher sampling rate in order to obtain the data samples needed to form the moment slices 
in accordance with (2.30) while satisfying (2.28) and (2.29). 
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V. CONCLUSIONS 


A. DISCUSSION OF RESULTS 

The 1-1/2 Dips method handles both signal and noise power differently than the 
spectrogram. With respect to signal power, the maximum processed signal power is the 
same for either method. But, because 1-1/2 Djps signal power is modulated, the average 
signal power is either a third of the spectrogram’s if the signal is real valued, or half the 
spectrogram’s if the signal is complex valued. With respect to noise power, the methods 
exhibit estimate variances that depend differently on the size of the data window and the 
input variance. Whereas the spectrogram’s variance is a function of input variance 
squared and window length squared, 1-1/2 D^pg’s variance depends upon input variance 
cubed and window length. Consequently, 1-1/2 Djpg is metre resistant to input noise than 
the spectro^am if the noise level is low enough and the window length is large 
enough. Conversely, the spectrogram handles noise better if either the window is too 
small or the input noise is too great. When die effects of signal and noise arc combined, 
the methods are compared by determining how many noise standard deviations the average 
signal power is above the mean noise level. This measure represents the method’s 
processing gain. 1-1/2 Dips ® larger processing gain than the spectrogram at high SNR 
for any reasonably sized data window. As SNR decreases, the processing gain advantage 
enjoyed by 1-1/2 Dips diminishes. At a sufficiently low SNR value that depends on the 
size of the data window, the spectrogram becomes the better estimate. As window length 
increases, 1-1/2 Djpj’s processing gain improves relative to the spechogram’s. 

As implemented in this study, the spectrogram and 1-1/2 Dips ®re superior to the 
bispectrum at detecting multi-component stationary signals in white Gaussian 
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noise. However, the frequency bin averaging scheme used in the simulations is not as 
favorable to the bispectrum due to the nature of its frequency-frequency representation. 
The bispecfrum does have an advantage though over the spectrogram and 1-1/2 Djpj in its 
ability to distinguish between unrelated signals, harmonic signals, and quadratically phase 
coupled signals. The only limitation is that quadratic phase coupling is not detectable if 
the frequencies of the coupled signals are harmonics of the same fundamental frequency. 

IHOMS usefulness is limited to detection of linear chirps having different slew 
rates. This method is considerably more costly than the spectrogram and 1-1/2 D^ps in 
terms of computations, especially for noisy signals. 

B. SUGGESTIONS FOR FUTURE STUDY 

Future work on two aspects of 1-1/2 Dip^ might prove worthwhile. The first is 
finding an optimum window and/or smoothing function. Rectangular windows, 
Hamming windows, and boxcar smoothing have been used in ONE_HALF to date. Since 
1-1/2 Dipj is a degenerate form of the bispectrum, a window function satisfying bispectral 
window requirements might prove to be better window. The second area warranting 
further study involves the inherent modulating effect of the 1-1/2 Djpg algorithm. The 
modulation contains useful signal information if it could be extracted from the 1-1/2 Djpj 
representation. 

The logical next step for future work on the bispectrum is to extend it to a time- 
frequency representation. Implementing an alternate bispectral calculation method 
relying on polar rasters looks like a promising starting point. Evaluation of parametric 
methods and comparison to the non-parametric methods of computing the bispectrum 
might prove worthwhile. In addition, better averaging schemes for the bispectrum 
frequency-frequency surface could be developed. 
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APPENDIX A: VARIANCE OF THE 1-1/2 Djps ESTIMATE 


The 1-1/2 Djpj spectrum is computed using (2.24). The variance of this estimate is 

derived here for a zero mean i.i.d. Gaussian input, (0, a^). Results are obtained for 

both real and complex valued noise. A rectangular window, i.e. w(^)=l for 0 ^ ^ V - 1, 
is assumed for these calculations. 

The mean of the 1-1/2 Djpj estimate is; 

N-l 

F {1-1/2 Dips) - 2 X ^ ^ in)\^x* in-k) + \x (n) (n + k) } exp (-joik). 

^k-0 

Since the expectation operation is distributive. 


1^-1 

£{1-1/2 Dip,]= 2 Y, («-^)] +E[\x{n)\'^x{/i + k)]} expi-jcak) . 

^k-0 

Let terml = E { |a: (n) |^a:* (n - /:) } and term2 = E {\x in)\'^x (n + k)}. Expand¬ 
ing term 1 yields 

term! = E {x’^ (n)x(n)x* {n - k)} . 


For k=^0, 

terml = E {x^ (n)xin)} E {x* (n-k)} 

= G^Ll 
= 0 * 

since the input samples are independent and the mean of the input sequence, |X^ = 0. 


For it =5 0, 

terml = E {x* {n)xin)x* (n)}. 
Assuming that the real and imaginary parts of x^n) are independent. 
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jc^(rt),x.(/i) ~N(0, ^aj), 

then terml can be rewritten as: 

terml= E { [x^ (n) -jx^ (n) ] [x^ (n) +jXf (n) ] [x^ (n) -jx^ («)]}, 
which simplifies to, 

terml= {E [x^ (n) ] + £ [x^ (n) ] £ [x? (n) ]}-;{£ [x] in)]+E [x^ (n) ] £ («) ] } 

after cross-multiplying and regrouping the real and imaginary parts. The terms in the last 
expression involving the product of a first moment and a second moment are zero because 
the first moment is zero. The third moment quantities in the same expression are also 
zero because the processes are distributed normally about a zero-mean. The final result 
is that terml is zero for any value of k. The same is true for term2. Consequently, 

1^-1 

£ {1-1/2 Dip^} - 2 ^ {terml + terml] exp (-jcok) 

= 0. 

The second moment of the estimate is: 

£{|l-l/2 Ap/} = £{1-1/2 1-1/2 CD)}. (A.l) 

where, 

-(m) 1^“^ 

l-m Dips = 2 {\xin)\^x* (n-m)+\x(n)\h(n + m) ] exp{-j(iim), 

m-O 

and 

l-ip. Dips (n,( 0 ) = 2 H, • 

j ■ 0 

After multiplying the last two equations together, collecting terms, and distributing the 
expectation operator, (A. 1) becomes: 
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I 2 J i * 

£{ 11 - 1/2 ip ~ ^ ^ {terml + term!+ term3 + term4} exp (‘~'j(a(m-s)) 

ntmOSmQ 

(A.2) 

where terml = E [\x{n)\^x* (n~^m)x{n-s)], 

term! = E[\x(n)\^x* (n-m)y^ (n + s)] 
term?) = E[|A:(/i)|^x(n + m)jc(n- 5 )], 
and term4 = E{\x{n)\^x{n + m):)t {n-¥ s)]. 

The first element in each of the terms above represents a product of four elements: 

\x(n)\^ = X* (n)x(n)x* (n)x(n). 

Terml, through termA are considered separately for all possible m and s index 
values. Two different mixed-moment relationships for a zero-mean Gaussian process are 
needed. For a real, zero mean Gaussian process, the following relation applies [Ref. 7]: 

E » E {X 2 X 2 } E {X 3 X 4 } +£ E {X 2 X 4 } (A.3) 

-j-£ 1 x^X 4 } £ {X 2 X 3 } . 

For a complex zero mean Gaussian process, a slightly different relation holds [Ref. 11]; 

E 1 X 2 X'*' 3 X 4 } - E {x‘^iX 2 } E {x* 3 X 4 } +£ {x“^lX 4 } E {x’^ 3 X 2 } . (A.4) 

Terml is evaluated first: 

1. For m 0 and s 5 ^ 0 , 

terml = £ [x'*'(n)x (/j)x* (n)x(n) ] £ [x* (n - 7 n)x (n-s) ]. (A.5) 

Applying (A.3) to a real process, rer/nl reduces to: 
rerml^ = 3 [£ {x(n)x(rt)} ]^£ {x(n-m)x(/j- 5 ) } (sincei.i.d.) 

9 2 

= (since zero-mean) 
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= 3a^ ■ a^6 (5 — m) (since i.i.d.) 

= 3a^h(s-m). 


(A.6) 


For a complex process, (A.4) is employed on (A.5) to achieve a similar result: 

terml^ - 2 [E{xin)x* (n)}]^E{x* {n-m)x{n-s)} 

- 2iol)^R^im-s) 

= 2a^ ■ a^6 (m-s) 

= 2op(m-s). (A.7) 

2. If m = 0 and s*0, 

terml^ = E {x* (n)x(n)x* {n)x{n)x* (n)} E {xin-s)} (A.8) 

= 0 (zero-mean process ^ £ {x (n - 5 ) } = 0) 

= terml^. 

This result is the same for both a real process and a complex process. 

3. If m 0 and 5 5 = 0, 

terml^ = E {x* (n)x(«)x* (/i)x(/i)x(n) } E {x* (n-m)} (A.9) 

= 0 (zero-mean process => E {x* (n-m)] =0) 

= terml^. 

4. If m = 5 = 0, 

terml^ ~ E {x* {n)x(n)x* in)xin)x* (n)x(n)} . (A.IO) 

Assuming that the real and imaginary parts of the process are independent 
implies that, 

x^(n),x,.(n) ~Ar(0, (A.ll) 

Dropping the index notation in (A.IO), breaking up each term into its real and 
imaginary parts, and grouping like terms produces: 

/crml^ = £{ (x^-7x.)^(x^-l-yx.)^} , 
which, after multiplication and collection of terms, becomes; 
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terml^ = +3£{a;^}£{4} +3E {x^^} E {x^} +£{xf}. (A.12) 

The second order moments in (A.12) are simply the variance of the real and 
imaginary parts of the process. The fourth order moments are computed using 
(A.3). The sixth order moments are obtained by computing the sixth deriva¬ 
tive of the moment generating function for a zero-mean Normal process [Ref. 
17]: 


E{p°} s= (0 = 

= 15o«; 


p 


where the subscript p denotes either or x, since their variances are equal. For 

a complex process, in terms of this variance, (A.12) becomes: 


termf = 15a6 + 3(ap (Sa^) +3(3a4) (ap + ISaJ 

= 48a®. 

p 

From (A. 11), = i:o\. Using this relation, (A. 13) becomes: 

P 2, * 


termv 


ra2f 

= 6a®. 


(A. 13) 


(A. 14) 


If the process is real, the x,- terms in (A.12) are zero and x^ terms become x: 


terml^ — £ {x®} 


Evaluating the sixth moment using the moment generating function yields: 


terml = 15a®. 


(A.15) 


Term2 was similarly analyzed. For the m ^ 0 and s=^0 index combination, the 
second moment term in the equivalent (A.6) and (A.7) is: 

£{x* (n-m)x(n + s)} - a^6(-m-s) 

= 0 , 
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because the delta function is always zero for non-negative indices, hike terml, term! is 
zero whenever just one index is zero because the process is zero mean. If the process is 
i: real and both indices are zero then term2 becomes a sixth order moment: 

terml^ ^ E {x (n) X (n)X (n) X (n) X (n) X (n) } 

= 15a^ 

If the process is complex, four terms in the expectation are conjugated while two are not. 
The result is a final expected value of zero since the process is zero mean Gaussian with 
independent real and imaginary parts: 

I terml^ = E {x* (n)x(n)x* {n)x(n)x* (n)x* (n)} 

= E{ix^-jxi)^ix^+jxi)^} 

~ 0 . 

Consequently, term! is zero for aU index combinations if the process is complex and is 
non-zero only for the 77Z = s = 0 condition if the process is real. 

I Evaluation of termS produces the same result obtained as for ternO., and evaluation of 

termA produces the same results obtained as for terml. Table A.l, and Table A.2 
\ . summarize term evaluation results for the real and complex processes, respectively. 


TABLE A.1: SUMMARY OF TERMS FOR REAL PROCESSES 


index 

condition 

terml 

term2 

terms 

term4 

m^O and 
^=5^0 

3a^b (s - m) 

0 

0 

3a^& (5 - m) 

m = 0 and 

0 

0 

0 

0 

m^O and 

5 = 0 

0 

0 

0 

0 

o 

II 

cn 

II 

s 


15a« 

15a‘ 

15 a® 

X 
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TABLE A.2: SUMMARY OF TERMS FOR COMPLEX PROCESSES 


index 

condition 

terml 

term2 

terms 

term4 

m*0 and 
5^0 

2a®6(5-m) 

0 

0 

2a®6(5-m) 

m = 0 and 
s^O 

0 

0 

0 

0 

m*0 and 

5 = 0 

0 

0 

0 

0 

m = s = 0 

6^' 

0 

0 

6af 


To obtain the second moment of the 1-1/2 Djpj estimate for a real process, the terms 

shown in Table A.l are substituted into (A.2). The m = s ^ 0 indices are taken out of 
the summation to obtain; 


1 


N~1N~1 


{\lA/2 Dip,\ }itEAL= 2(3ap6(5-m)exp(-ytD(m-5))} 

m - li - 1 


1 


N-1 


= ^ {60a® + ^ 2 (30p exp {-jcam) } (sifting theorem) 


m- 1 


= ^{60aJ+ (iV-l)6ap 

54 + 6A, 

= --^(ap; 


(A.16) 


where © = 0 was assumed going from the second line to the third line in order to account 
for the maximum variance contribution from the sum term. 

Following the same steps for a complex signal, (A.2) becomes; 

E {\l-iP. dip^\ } COMPLEX ~ • (A^17) 
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The variance of the 1-1/2 Dzpj'estimate is given by; 

var {1-1/2 = £{|l-l/2£>J^ - [£{1-1/2 

Note that the variance is simply equal to the second moment since the mean of the 
estimate is zero. Therefore, the variance of the 1-1/2 estimate is given by (A. 16) for a 
real zero-mean independent Gaussian process. If the process is complex zero-mean 
independent Gaussian with independent real and imaginary parts, the estimate’s variance 
is given by (A. 17). 


s 


APPENDIX B: THEORETICAL 1-1/2 Djps SIGNAL POWER 

This Appendix first shows how the theoretical 1-1/2 Djps spectrum is calculated for an 
input signal consisting of a single real valued sinusoid. The computation is then repeated 
for a complex valued sinusoid. Both calculations are made using a rectangular window in 
(2.24). 

A. REAL VALUED SINUSOID CALCULATION 
The real valued sinusoidal input signal is 

x(n) = cos (0n), (B.l) 

where, 


0 = 2n~, (B.2) 

is the digital frequency of the signal. Substituting (B.l) into (2.24) using a rectangular 
window yields: 

v.l 

1-1/2 Djpj (n, ®) = 2 “ ^') ] + cos [0 (n + ^) ] } exp (-jwk ). 

kmO 

(B.3) 


After applying the trigonometric cosine product identity, 

cos (y4-B) + cos (A+B) ~ IcosAcosB, (B.4) 

and pulling the cosine cubed term out of the summation since it is not dependent on k, 
(B.3) becomes; 

iV-l 

1-1/2 D.pj,(n,a)) ~ ccs^ (dn) ^ cos (dk) exp (~Jcok). (B.5) 

it - 0 

The summation portion of (B.5) is the Fourier transform of the sinusoidal signal. 
Evaluating the transform at co = 0 and expressing 0 as shown in (B.2), the Fourier 
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transform can be expressed in terms of delta functions; 


N-l 


^ COS {2n^k) exp i-j2n~k) = ^ { 6 (®-m)+ 6 [c 0 “‘(A^-m)] } . 
kmO 

With (B.2) and (B.6), (B.5) now becomes: 

W 3 yn 

1-1/2 D.pjjCn,©) = ^cos (27t^n) {6(to-m)+6[®-(//-m)]} . 




N 


(B.6) 


(B.7) 


B. COMPLEX VALUED SINUSOID CALCULATION 
For a complex valued sinusoid input, 


m 


(B.8) 


x(n) 5= exp(J2n-n) 

= expijdn), 

the 1-1/2 Djpj estimate calculated by (2,24) using a rectangular window is; 

1^-1 

1-1/2 Djp^ in, ©) = - \exp ijQn)\^ [exp [-;0 (n -1:) ] + exp [jd (n + k)]} exp {-jcok) 

(B.9) 

After factoring out the exponential terms that are not a function of k, and recognizing that 
the magnitude squared term is simply equal to one, (B,9) reduces to: 


1 


N-l 


1-1/2 Djpj (n, ®) = 2 { exp [/6n] + exp [-/On] } exp (jdk) exp (-jcak). (B.IO) 

i-o 


With the Euler identity. 


1 


2 L/®«] + exp [-jOn] } = cos (0n), 

and (B.2), (B.IO) becomes: 

N-l 

1-1/2 Djpj in, ®) = cos i2n~n) ^ exp iJ2n~k) exp i-ja)k ). 


(B.ll) 


(B.12) 


ifc-O 
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The summation term in (B. 12) is the Fourier transform of a complex sinusoid which is 
expressible in terms of a delta function yielding: 






APPENDIX C: COMPUTER PROGRAMS 


The four extrinsic Matlab functions contained in this Appendix are written for Matlab 
version 4.0. If the functions are used with earlier versions of Matlab, contour plots and 
mesh plots will not be oriented correctly. 



%function [P,fTeqindex,tiniemdex] = one_half(data,wintype,winlen,step); 

%This function will calculate the 1 1/2 Dips spectral surface. The 1 1/2 Dips 
%surface characteristics are determined by the selection of window 
%type (wintype), window length (winlen) and the distance that the window 
%is moved through the data sequence (step). The magnitude of the positive 
%half of the 1 1/2 D spectral plane is returned in the “P” matrix. 

%Outputs “timeindex” and “freqindex” aid in axis labeling. 

%The inputs are: 

%data - The input observations vector. The length should be a power of 2. 
%wintype; ‘0’ Rectangular Window 
% ‘ r Hamming Window 

%winlen - The desired width of the window, normally half of the input 
% length. 

%step - Time step desired, can be ‘1’ or a multiple of ‘2’ 

%prepared by Karen A. Hagerman, 06 May 1992. 

%modified by Jeff McAloon, 01 June 1993. 

function [P,freqindex,timeindex,svect]=one_haIf(data,wintype,winlen,step) 

[datarows,datacolumns]=size(data); 
if datarows~=l 
data=data.’; 
end 

siglen=length(data); 

96 







if wintype==0 

win=ones(winlen-1,1); 
elseif wintype==l 

win=hamming(winlen-1); 
end 

W=[win(winlen/2:-1:1)]; 
x-[zeros(l,winlen) data zeros(l,winlen)],’; 
p=zeros(siglen/step,winlen); 

mdex=l; 

for n=winlen+l :step:siglen+winlen'Step+l 
xt=x; 

ll=max((winlen+l),(n-(winlen/2-1))); 
ul=min((n+(winlen/2-1 )),(siglen+winlen)); 
ctp=xt(ll;ul); 
mtp=mean(ctp); 
xt(ll;ul)=ctp-mtp; 

xdt=xt((n'(winlen/2-1)); (n+(winlen/2-1))); 
Xn=[abs(xt(n))^2;abs(xt(n))^2]; 

Xin=[conj (xt(n:-1 :n-(winlen/2-1))) xt(n:n+(wmlen/2-1))]; 
product=((Xm*Xn).*W).’; 
product=[product 0 conj(product(winlen/2:-l:2))]; 
pt=0.5*fftshift(fft(product)); 
p(mdex,:)=pt; 
index=index+l; 
end 

P=abs(p(; ,winlen/2+1: winlen)); 
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[prow,pcolumn]=size(P); 
freqindex=tO:pcolumn-1 ]; 
tiineindex=[ 1 .prow]; 
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% function [P.ffeqaxis.timeaxis] = spectro(data,wintype,winlen); 

% 

% This function calculates spectrogram of the supplied sequence, 

% “data”. The user must specify: 

% “wintype” - “0” for a rectangular window. 

% “1” for a Hamming window. 

% “winlen” - The desired length of the data window, 

% The time step is fixed at one and the spectrogram is calculated with 
% non-normalized periodograms. The time-frequency surface is returned 
% in the “P” matrix. The columns of “P” are the frequency bins while the 
% rows are the time steps. Time-frequency axis labeling vectors, 

% “freqaxis” and “timeaxis” are also returned to aid in the plotting 
% of results. 

% prepared by Jeff McAloon, 01 June 1993. 

function [pow,freqaxis,timeaxis] = spectro(data,wintype,winlen) 
x=data; 

xlen=length(x); 

N=^inlen; 

[row, col]=size(x); 
ifrow>l 
x=x.’; 
end 

x=[2eros(l,N/2) X zeros(l,N/2)]; 
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if wintype == 1 
win=haniiningCN) ’; 
elseif wintype == 0 
win=ones(l,N); 
end 

for ind=l:xlen 

xseg=x(ind; (N-1+ind)); 
xseg=xseg-mean(xseg); 
SP=(abs(fft(xseg.*win,N))) 
POW(ind,:)=fftshift(SP); 
end 

pow=POW(:,N/2+l;N); 
timeaxis=0: xlen-1; 
ffeqaxis=0;N/2-l; 


% function [bis, baxis] = indbis(data,datseg,lag,tlen,wm) 

% This function computes the bispectrum of real and complex valued data 
% sequences by the indirect method. The tri-correlation sequence is 
% computed for a specified number of segments. These segment sequences 
% are then averaged to obtain the final tri-correlation sequence. The 
% three available tri-correlation window options are the unit hexagonal 
% window, Parzen’s window, and the optimum window (also known as Sazaki’s 
% minimum bias window). The required function call arguments are; 

% “data”-> input data sequence vector. 

% “datseg”-> number of samples to be used in each segment. 

% “lag” -> number of lags to be used in the computation of the 
% tri-correlation sequence. 

% “tlen” -> square dimension of the two dimensional FFT to be used 
% on the fmal form of the tri-correlation sequence. 

% “win” -> “0” for unit hexagonal window 
% “1” for optimum window 

% “2” for Parzen’s window. 

% The function returns the complex bispectrum in the “tlen-by-tlen” sized 
% matrix “bis”. A vector “baxis” is returned that can be used to label 
% both bispectrum frequency axes. 

% prepared by Jeff McAloon, 01 Jun 1993. 

function [bis, baxis] = indbis(data,datseg,lag,tlen,win) 

% Initialize parameters and reorient data if necessary: 
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N=length(data); 

Ns=fix(N/datseg); 

[rd, cd]=size(data); 
if cd> 1 
x=data.’; 
else 
x=data; 
end 

x«[0;x-mean(x)]; 

Rt=zeros(2*lag+1); 
rx=zeros(2*lag+1); 

% Compute tri-correlation sequence: 
for i!=l:Ns 
for j=l:datseg 

y(j)=x(j+(i-1 )*datseg); 
end 

y=y-mean(y); 
for m=l;(2*lag+l) 
for n5=l:(2*lag-f-l) 

sl=max([l,-(m-(lag+l)),-(n-(lag+l))]); 
s2=min([lag,lag-(m-(lag+1 )),lag-(n-(lag+1))]); 
n=0; 

fork=sl:s2 

^=^+conj(y(k-^l))*y(k^-l+(m-(iag+l)))*y(k-^l-^(n-(lag+l))) 

end 

rx(n,m)=T; 

end 
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end 


Rt=Rt+rx; 

rx=[]; 

end 

R=flipud(Rt/Ns); 

R=Rt/Ns; 

% Determine tri-correlation window function: 

if win = 0 % unit hexagonal window 

W^ones(size(R)); 

elseifwin==l % optimum window 

itm=0; 

for m=-lag:lag 
itm=itm+l; 

wm=abs(sin(pi*m/lag))/pi + (l-abs(m/lag))*cos(pi*m/lag); 
itn=0; 

for n=-lag:lag 
itn=itn+l; 

wn=abs(sin(pi*n/lag))/pi + (l-abs(n/lag))*cos(pi*n/lag); 
wmn=abs(sin(pi*(m-n)/lag))/pi + (l-abs((m-n)/lag))*cos(pi*(m-n)/lag); 
W(itm,itn)=wm*wn*wmn; 
end 
end 

elseifwin==2 % Parzen’s window 

itm=K); 

for m=-lag:lag 
itm=itm+l; 
if abs(m) <=: lag/2 
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wm=1-6* (abs(m)/lag)^2+6*(abs(m)/lag)''3; 

else 

wm=2*(l-abs(m)/lag)^3; 
end 
itn=0; 

for n=-lag:lag 
itn=itn+l; 
if abs(n) <= lag/2 

wn=l-6*(abs(n)/lag)'^2+6*(abs(n)/lag)^3; 

wmn=l-6*(abs(m-n)/lag)'^2+6*(abs(m-n)/lag)^3; 

else 

wn=2* (1-abs (n)Aag)''3; 
wmn=2*(l-abs(m-n)/lag)^3; 

end 

W(itm,itn)=wm*wn*wmn; 

end 

end 

end 

% Compute bispectrum: 
bis=fftshift(fft2(R*W,tlen,tlen)); 
baxis=(-tlen/2): (tlen/2-1); 
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% function Amtrx = gen_a(al_elem); 


% This function computes the parameter vectors needed to compute a fourth 
% order moment slice in the IHOMS time-frequency method implemented by 
% the extrinsic matlab function ATH_IMS.M. The first element, al, of 
% each IHOMS parameter vector is specified in the function call as the 
% vector, “al„elem”. The length of “al_elem” is the number of distinct 
% parameter vector sets (i.e.one set would be [1 al a2 a3]) returned as 
% columns in the matrix “Amtrx”. This also equals the number of time- 
% frequency surfaces summed by IHOMS to form the final time-frequency 
% representation. 

% prepared by Jeff McAloon 01 June 1993. 

function Amtrx = gen_a(al_elem) 

% Verify input vector a column vector: 

[rin, cin]=size(al_elem); 
if cin > 1 

alv=al_elem.’; 

end 

% Iteratively solve for parameter vectors: 

A=0; 

for ps=l:length(alv) 
a3v=-5:l:5; 

a2a=100; a2b=l; i=l; md=99; 
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while (abs(a2b-a2a) >= 0.1) & (i <= length(a3v)) 
a2a=l-alv(p)-a3v(i); 
a2b=real(sqrt(l-alv(p)'^2-a3v(i)^2)); 
if abs(a2a-a2b) <= md 
indmin=i; 
end 

md=mm(ind,abs(a2a-a2b)); 

id=max(2,indinin); 

i=d+l; 

end 

a3v=a3v(id-l):0.0001:a3v(id); 
a2a=l; a2b=2; i=0; 

while (abs(a2b-a2a) >= 0.0001) & (i <= (length(a3v)-l)) 
i=i+l; 

a2a=l-alv(p)-a3v(i); 

a2b=sqrt(l-alv(p)''2-a3v(i)''2); 

end 

al=alv(p); a2=a2a; a3=a3v(i); 
a=[lala2a3]; 

A=[A,a.’]; 


Anitrx=A; 
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% function surf = ath_ims(M,u,b,N,SNR,A.g) 


% This function computes the IHOMS time-frequency representation of 
% multi-component chirp signals in white Gaussian noise. Input arguments 
% are: 

% M -> number of chirp components. 

% u -> vector of chirp starting frequencies. 

% b-> vector of halved chirp slew rates. 

% N -> data record length. 

% SNR -> desired SNR. If noise free signal desired enter “99”. 

% A -> matrix whose columns are the “a” parameter vectors needed to 
% form the moment slice. See extrinsic Matlab function GEN_A 
% generate this matrix. 

% g-> kernel function parameter. 

% The moment slice is computed using 100 lags. The IHOMS surface is 
% plotted along with a 1-D radial maxima plot that aids in locating and 
% characterizing the chirps present in the signal. In addition, the 
% function returns a 100 X 100 matrix, “surf”. The columns are time/lag 
% and the rows are frequency bins. 

% prepared by Jeff MeAloon, 01 June 1993. 

function surf - ath_ims(M,u,b,N,SNR,A,g) 

% Initialize parameters: 
randn(‘seed’,0) 

if SNR ==99 % noise free 
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GG=0; 

else 

GG=1; 

end 

Z=sqrt(10^(GG*SNR/10)); 

inaxlag=100; 

j=sqrt(-l); 

[rA,cA]=size(A); 

SURF=zeros(100); 

% IHOMS algorithm: 

fork=l:cA 

a=A(:,k); 

kn=k 

TF=D; 

for n=0:maxlag-l 
for lag=0:maxlag-l 

no=GG*randn( 1 ,length(a)); 

term(l)=0; 

forp=l;M 

term(l)=term(l)+Z*exp(j*2*pi*(((lag)/N)*u(p)+(((lag)/IS0''2)*b(p))); 

end 

term( 1 )=term( 1 )+no( 1); 
for i=2:length(a) 
arg=n+a(i)*lag; 
if arg < 0 
term(i)=0; 
else 
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tenn(i)=0; 
for q=l:M 

term(i)=term(i)+Z*exp(j*2*pi*(((arg)/N)*u(q)+(((arg)/N)^2)*b(q))): 

end 

term(i)=term(i)+no(i); 

end 

end 

kern=exp(-g*lag^2); 

ni(lag+l)=kem*conj'(tenn(l))*tentt(2)*term(3)*term(4); 

end 

TF(n+l.:)=fft(m); 

end 

SURF=SURF+TF; 

end 

% Orient matrix for MATLAB 4.0 contour plot: 
surf=rot90(fUpud(abs(SlJRF)),-1); 

% Plot IHOMS surface: 

t=0:maxlag-l; 

f=0:maxlag-l; 

subplot(121),contour(t,f,surf,2),title(‘IHOMS surface’),... 
xlabel(‘tiineAag’),ylabel(‘frequency’) 

% IHOMS 1_D search routine and plot: 

theta=0:0.01:pi/2; 

for in=l:length(theta) 

D(in)=0; 
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forr=30:90 


t=l+r*cos(theta(m)); 

f=l+r*siii(theta(in)); 

D(m)=D(in)+surf(f,t); 

end 

end 

subplot(122),plot(theta,D),xlabel(‘theta in radians’), 
ylabel(‘D(theta)’),title(‘l-D radial maxima plot’) 
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